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Creep Collapse of Long Cylindrical Shells 
Under High Temperature and External 
Pressure’ 


THEIN WAH* anp R. KIRK GREGOR Y** 
Southwest Research Institute 


Summary 


This paper describes theory and tests of the creep collapse 
of long thin aluminum-alloy cylinders under external radial 
pressure. Steady-state creep is assumed in the theoretical deriv- 
ation. The test temperatures were between 300° and 500°F. 
The collapse time for each cylinder was calculated theoretically. 
Agreement between theoretical and test results was fair. 


Symbols 


mean radius of cylinder (in.) 

twice the thickness of membrane facings of sandwich 
(in.) 

natural base of logarithms 

nondimensional deflection 

initial nondimensional deflection at ¢ = O after load 
application 

initial nondimensional deflection before load applica- 
tion 

thickness of cylinder wall (in.) 

thickness of solid wall of cylinder (in. ) 

thickness of sandwich wall (in. ) 

applied pressure (psi) 

elastic buckling pressure, (psi) 

gas constant, taken as 2 cal./mole/°K. 

constant, per hour per °K. 
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temperature (°K.) 

time (hr. ) 

collapse time (hr. ) 

deflection (in.) 

initial deflection after load application (in.) 
initial deflection before load application (in.) 
Ger/Q 

(&/o0)foo 

activation energy (cal./mole) 

(¢/a0)(B/(B — 1)] 

strain 

curvature 

Poisson’s ratio 

mean stress (psi) 

constant in creep law (psi) 

cylindrical coordinates 

temperature constant defined by Eq. (18) 


(1) Introduction 


Sus ADVENT of certain high-temperature operating 
conditions has rendered urgent studies on the be- 
havior of various structural configurations under 
environmental conditions not often encountered on 
earth. For example, a body moving rapidly through 
air is subjected to aerodynamic heating, which causes 
the temperature of the body to rise. The increase in 
temperature affects the material of the body and, in 
particular, causes it to creep—that is, to continuously 
change its dimensions under load. At elevated tem- 
peratures, all known materials experience a certain 
amount of creep, and this must be taken into account 
in strength calculations. 

In the case of structures which may become un- 
stable under load, the creep of these structures at 
elevated temperatures may be sufficient to cause them 
to deform excessively. When this occurs, the struc- 
ture is no longer able to fulfill its intended function 
and it has reached the end of its useful ‘“‘lifetime.”’ 
This is usually referred to as the critical time, ¢.,. 
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Fic. 1. Schematic of sandwich wall. 


Most of the studies in this area have heretofore been 
directed at predicting the effect of creep on the life of 
columns and plates under compression. The “‘life- 
time’ of a cylindrical shell under external radial pres- 
sure and high temperature has received relatively little 
attention and the authors have been able to uncover 
only four papers in English (references 1, 2, 3, and 4) 
dealing specifically with this problem. 

A long cylinder, subjected to a constant external 
radial pressure q less than the static buckling pressure 
Jer, and to a reasonably high temperature for a suffi- 
cient length of time, will eventually collapse, however 
small the constant pressure g. Analysis shows that 
the mechanism of this collapse is dependent on the 
presence of initial deviations from perfect circularity 
in the shape of the cylinder cross section. Conse- 
quently there exists, from the very start, a differential 
stressing of the fibers of the cylinder wall and there- 
fore a differential creep rate which eventually causes 
collapse. It follows from this that the creep-collapse 
problem cannot legitimately be classed a characteristic- 
value problem in the mathematical sense. The struc- 
ture fails as a result of continually increasing deforma- 
tions, and not because at the so-called ‘‘critical time’’ 
a bifurcation of equilibrium positions becomes possible. 
The terms ‘‘creep buckling’? and ‘creep instability” 
often encountered in the literature must therefore be 
interpreted in this restricted sense. 


(2) Scope of the Present Investigation 


This study is limited to the investigation of the creep 
collapse of a long thin tube with a known out-of- 
roundness, which is subjected to an external radial 
pressure g at a constant temperature between 300° and 
500°F. In particular, an attempt is made to predict 
theoretically the ‘“‘collapse time’ or ‘“‘critical time” 
t-, of such a cylinder under given conditions of temper- 


ature and pressure. For analytical purposes, then, it 
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is assumed that the cylinder is infinitely long. The 
theoretical derivation is given in Section 3. 

A number of experiments were conducted to verify 
the theory, and the test results are summarized and 


discussed in Sections 4 and 5. 
(3) Theoretical Analysis 


The analysis of the creep collapse of long tubes re- 
quires certain simplifying assumptions at the outset. 
Cnce these are made, the analysis can proceed in a 
manner very similar to that used in deriving the critical 
pressure of a long cylinder in the well-established case 
of elastic buckling. 

These assumptions are as follows: 

(1) The solid wall of the tube can be replaced by a 
sandwich wall as suggested by Hoff, Jahsman, and 
Nachbar! and by Sundstrom.” The sandwich wall is 
supposed to be such that the outer and inner sheets 
can resist only membrane stresses, whereas the core 
material takes no normal stresses and has infinite 
resistance to shear deformation. It is realized that 
this type of wall represents the solid wall only approxi- 
mately, but for thin walls the approximation should be 
close. The sandwich-wall tube is illustrated in Fig. 1. 

(2) While the tube is probably in a state of plane 
strain, it is sufficient to assume a state of plane stress. 
The principal reason for this assumption is the diffi- 
culty of adapting the creep law used for a state of bi- 
axial stress. A supplementary analysis shows, how- 
ever, that an ap) roximate extension of the creep law 
to apply to biaxial stress results only in a minor modi- 
fication of the resulting equations. 

In Fig. 2(a), the coordinate system used is shown, and 
This condition will 


plane strain means that «, = 0. 
be most nearly true in the center portions of the cylin- 
der, but will be violated near the ends. Plane stress 
means that o, = 0, a condition that is true only at the 
ends of the cylinder, where the support condition is 
such that no axial stress exists. 

Consider the cross section of an initially out-of- 


round cylinder, as shown in Fig. 2(b). This out-ol- 














Fic, 2(a). Coordinate system. 
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Fic, 2(b). 
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roundness makes the tube slightly elliptical and is the 
form of initial deformation most commonly encoun- 
tered. 

The initial deviation from circularity is then defined 


by the expression 
Wo = Hf cos 20 (1) 


where // is the thickness of the wall. The symbol h, 
will be used for the thickness of the sandwich wall 
(Fig. 1) and f# for the equivalent solid wall. When an 
expression applies to both solid and sandwich walls, 
the symbol /7 will be used, as in Eq. (1), where foo is 
the radial deviation from circularity at 6 = 0, 7, non- 
dimensionalized by the thickness /7/. 

Suppose next that a uniform external pressure gq, 
less than the elastic critical pressure q,,, is applied to 
the cylinder. On page 223 of reference 5 it is shown 
that the additional deflection as a result of the appli- 
cation of this pressure for a solid-wall tube is 


w’ = [q/(der — G) hf cos 26 (2) 


For a sandwich type of construction, the expression 
(2) will also hold if the value of g., is changed appro- 
priately. One may therefore write the total initial 


deflection at time ¢t = Oas 
Wii-0 = Wo = Wo + w’ = Hf cos 20 + 
[a/(Ger — 7) [foo cos 20 
or, 
wy = Hf, cos 20 (3) 


in which 


fo = foo/(1 — (9/ger) ] 
E a ' 
Ger * — — for the solid-wall tube 
al — 9°) @ 
3 E Ae... ; 
- — for the sandwich-wall tube 
4(1 — »*) a° 
The expression for g., in the case of the sandwich- 
wall tube may be deduced from the corresponding ex- 
pression for the solid-wall tube by noting that the 
centroidal moment of inertia per unit length of the tube 
is h?/12, whereas that of the sandwich tube is ,°b/4. In 
Eqs. (4), a is the mean radius of the tube cross section, 
E is Young’s Modulus, and »v is Poisson’s ratio for the 
material. , 
It is further shown in reference 5 that the deformed 
cross section of the tube may be regarded as the 
funicular curve for the bending moments in the wall as 


a result of the applied pressure g. One may thus write 


M, tu = Ma = qgawy 
or My, = gaHfy cos 20 (5) 
On the basis of Eqs. (3) and (5) it is now assumed that 


the radial deflections and bending moments at any time 
are given by the expressions 


w = Hfcos26 | 


” 
M, = qallf cos 264 (9) 
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in which f is now a function of time such that at ¢ = 0, 
f = fo. 

If w is positive inward, the curvature « of a ring is 
given by (see p. 206 of reference 5) 


xk = (1/a*)[w + (d°w/dé?) (7) 


Substituting the first of Eqs. (6) into (7), there 
follows 


k = —(3/a*)Hf cos 26 (8) 
Differentiating (8) with respect to time, one has 
k = —(3/a?)Hf cos 20 (9) 


in which the dot represents differentiation with respect 
to time and « stands for rate of change of curvature. 

It is noted that the assumed mode of deformation in 
Eqs. (6) implies a linear relation between moment and 
curvature (and moment rate and curvature rate) 
throughout the creep-collapse process. This is the 
natural result of the neglect of the superposed higher 
modes. Experimental data strongly suggest, how- 
ever, that the constraint imposed by this assumption 
is not a severe one while, on the other hand, it leads to a 
substantial simplification in the basic differential equa- 
tion (Eq. 16, below). 

The problem now reduces essentially to finding a 
relation between the moment \/, and the rate of 
change of curvature, «. This is done by writing an ex 
pression for the curvature rate in terms of the strain 
rates in the tube wall, and then making use of a creep 
law. These relations evidently depend on the type of 
construction and the creep law assumed. 

(3.1) Ideal Sandwich Construction-Plane Stress 

Referring to Fig. 1, the following equations of equi- 
librium may be written: 

a, = & + (2M,/bh,)\ 
a2 = & — (2M,/bh,)$ 


& = (1/2)(o, + o2) = ga/b (11) 


(10) 


\M, is considered positive when it produces a decrease 
in the curvature of the wall and compressive stress is 
taken as positive. The subscripts 1 and 2 refer to the 
outer and inner membrane, respectively. 

From the geometry of sandwich cylinder, Fig. 1, one 
niay write an expression for the curvature as follows: 


kK = —(e@ — &)/h, (12) 
Differentiating (12) with respect to time, one has 
kK = —(€& — &)/h, (13) 


While there are numerous phenomenological creep 
laws, it was decided that the formulations due to 
Elbridge Z. Stowell are theoretically the most accurate 
and have also been remarkably well verified by existing 
experimental data (references 6 and 7). For simplicity 
the steady-state creep law for the case of uniaxial 
stress, as given by reference 6, was chosen. This may 
be written in the form 
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Fic. 3. Collapse time for long tubes. 


1 1 {RT 
; ; (<.) + 2sTe~ S4/RY) sinh > (14) 


E 00 

where 

€ = strain 

E = elastic tensile modulus, psi 

o = Stress, psi 

7 = temperature, °K. 

s = constant, per hour per °K. 

AH = activation energy, cal./mole 

R= gas constant, taken as 2 cal./mole/°K. 

t = time 


In Eq. (14) the first term on the right-hand side is 
the elastic component of the strain rate due to a 
change in the stress. 

Writing Eq. (14) for the outer and inner fibers, one 
finds 


; d {a at whet. 9 
qo = ( ') + 2sTe~ S2/RY) sinh 


dt E a0 
(15) 


; d [a pas ci /RT) _- 02 
é& = - ( ) + QsTe~ GH/RY sinh — 
dt \E a0 


Substituting ‘or o; and o2 from Eqs. (10), and sub- 
tracting, there follows in view of Eq. (13) 


a-@ ._ 4M, “ 
Is poe Ebh,? 
dsTe~ GH/RD) 9 ; 


inh —~" cosh (16) 
h.K 8%: bh,oo a a 


Noting from the second of Eqs. (4) that the value of 
Ger for the sandwich-wall tube with v= 0 is given by 


Ger = (3/4) (Eh,*b/a*) 


and making use of Eqs. (6) and (8), one may write 
(16) in the following form: 


i a ) 
iia -l\d@ = 
3he?[(q/der) — 1]/ ‘ 
d(f cos 20) 


cosh (&/oo) sinh [2(¢/o0) f cos 26] -" 
where H has been replaced by /,. 
Integrating from ¢ = 0 to ¢ = ¢, one finds finally 
t cosh (¢/go) 
(ave 24/8") /2EsT)[(Ger/q) — 1] 
coth|[(a/a0)/o cos 26] 


og, . 
coth[(/oo)f cos 26] 


in which the initial condition f = fo when ¢ = 0 has 
been introduced. 
Using the notation 


r = 2.3026 ave” RY /2EsT\ 


18) 
log = logio f _ 


the above equation may be written 
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CREEP 


t 7 l 
t[(Ger/q) — 1] 


~ cosh (/a0) 
J a q 9 a P =a 
ylos coth — fy cos 20) — log| coth — f cos 2 


do 


The “‘critical time’ or collapse time will be defined 


as the time at which the deflection 

comes indefinitely large. There follows from 
ler log [coth(&/oo)fo| 

T((Ger/q) ~ By 


This is the basic equation proposed herei 


cosh(é/ ao) 


collapse time of an infinitely long tube. 
From Eggs. (19) and (20), one may write 


—- INI 


COLLAPSE OF 


f cos 20) maz be- 
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(ter — t) 
T((Ger/q) — 1] 


log [coth(¢/ao)/| 
' : (21) 


cosh( a/ oa ) 


Eq. (21) defines the time rate of change of the deflec- 


) (19) 


tion f and is suitable for experimental verification. 
(3.2) Solid Walls 


The extension of Eqs. (2()) and (21) to solid-wall tubes 


> (  * . ° . . or ~ 
Eq. (19) is immediate except for the interpretation of fo. From 
the equivalence of the moments of inertia and cross- 
(20) sectional areas of the sandwich- and solid-wall tubes, 


there follows: 


f he ; 
+ ae oe Teana Teoria = (bh,2/4) : (43/12) = 1:1 


Aas ‘SAolita = (b) : (h) <3 


TIAL OUT-OF-ROUNDNESS PARAMETER f, 
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Fic. 5. Error in ¢,, resulting from inaccuracy in foo (Or er). 


Thus 


h, = hvVf/3 


On the other hand, f (and fy) in Eqs. (20) and (21) is 
the out-of-roundness which has been nondimension- 


alized with respect toh,. Hence 


fo = wWo/h, = V3wo/h l ( 
s/o) = gas/hoy § 


bo 
bo 
Nahe 


The authors are indebted to an anonynious referee for 
clarifying some of the above relations. 

Figs. 3 and 4 show a plot of Eq. (20) with the left- 
hand side of (20) as ordinates and fj as abscissae. It 
may be remarked that since fj cos 26 is a constant for a 
fixed 6, the abscissae may equally well be taken as fo 
cos 20. These curves may be regarded as universal 
creep-collapse curves, as they are independent of tem- 
perature, 7 defining the temperature constant for a 
given material. Values of 7, for a given material, 
could be furnished in tabular form. 

It may be noted if the ordinates in Figs. 3 and 4 are 
made to give the left-hand side of Eq. (21), the abscis- 
sae will define f (or { cos 26)—1.e., the same curves may 
be read to yield the value of the deflection f at any 
time from ¢ = 0 to very near the theoretical collapse 


time f,,. 


(3.3) Ideal Sandwich Construction—Plane Strain 

An approximate adaptation of the creep law, Eq. 
(14), for the case of plane strain is obtained by writing 
itin the form 


Be 5 fea 
dt dt E | 


2sTe~ S4/RT) sinh Lao — (oe + 2)/3)} (23) 


a0 





The condition of plane strain may then be written 
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de, d Jo. = Vde\ 
= 7 


dts it z YT 
(og + o-)/3] 


r Oo: — 
sinh =() (24) 
00 


25 Te 7 (SH/RT) 


In Eqs. (23) and (24) the deviatoric stress takes the 
place of the stress in the second term on the right-hand 
side of Eq. (14). The first term on the right-hand 
side of Eqs. (23) and (24) follows from elasticity theory. 

On the physically reasonable assumption that the 
elastic and viscous components of the strains (and 
velocity strains) separately vanish, one may deduce 
from Eq. (24) that 


: = vog for the elastic term | 15) 
‘ ‘ F (20 
: = o,/2 for the viscous term} 
In view of (25), Eq. (23) may be written 
dé l— *\d 7 ,—(AH/RT) |: (09/2) 
= : (og) + 2sTe sinh 
dt E dt or 
(26) 


Comparing Eqs. (26) and (14) it is readily seen that the 
analysis for the case of plane strain, parallel to that for 
plane stress given in the previous section, would yield 
equations similar to (20) and (21), the only modifica- 
tion being that ¢ in Eqs. (20) and (21) would be replaced 
by «/2. It is noted, however, that Eq. (26) does not 
have a strong theoretical justification. 


(3.4) Sources of Error 
Apart from the errors introduced by the scatter of 
experimental creep data, which is considerable, pre- 
dictions of the critical time based on Eq. (20) may 
result in discrepancies because of the inaccuracies in the 
values of the elastic buckling pressure g.;, and in the 
measured value of the out-of-roundness parameter /oo. 
(These questions are discussed further in Section 5.) 
Error Introduced by q-;—In order to estimate the 
error introduced by an inaccuracy in the value of ¢-,, 
Eq. (20) is written in the form 
,AH/RI - 
log, qeoth : foo 
oi 


2EsT Siete |. 


where B = Ger/q 


Differentiating (27) with respect to 8, and dividing 
by ¢,,, one finds 


4, 278 
de ; ie esch - 
—-= — + —— — dg (28) 
log, | coth 
a= 1 
in which y = (&/o0)fo (29) 


From Eq. (28) it is seen that if an inaccuracy A@ oc- 
curs in the value of §, so that the true value of £ is 
(1 + A)8, then the corresponding change Af,, in f.; is 
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Fic. 8. Diagram of tape seal at end of specimen. 
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\C- Fic. 7. Details of test facility. Fic. 9. Mounted cylinder and post with deflectometers. 
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in which 


é ( —1 
y= 2 esch Pe | toe. (cotn 3 = :)| (31) 


Fig. 5 gives a plot of y against foo, for various values 
of the parameter 6 = (&/00)[8/(8 — 1)). 

Error Introduced by fy—In order to estimate the error 
introduced in ¢,, by an error in the measurement of 
foo, Eq. (20) is written in the form 


ter = A log,(coth foo) (32) 


in which A is aconstant and 


2GE)-tb-a 
oo \G — 1 eri) Wer 


Differentiating (32) with respect to foo and dividing by 
t.,, there follows 


(33) 


dtr 26 csch 26f00 
—_—_ =o —_— - , dfoo 
70 log,(coth 6fo0) 
1 Ate, 26(esch 26foo) 
or — see. rete (34) 
Afoo Fay log,(coth foo) 


In Eq. (34) At,, is the error in ¢,, as a result of an in- 
accuracy Af in foo. 

For convenience in plotting, Eq. (84) may be put in 
the form 


Abee/ ber = — dy(Afoo) (35) 
in which 
Y = 2(csch 25fo)/log,(coth 6 foo) (36) 


Fig. 5 gives a plot of foo against y, with 6 as a parameter. 
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(4) Creep Collapse Experiments 


(4.1) Testing Procedure 


Thirty tests were run on thin-walled tubes of alumi- 
num alloy 3003 H-14 to verify the theory. The tube 
was mounted on a steel shaft and supported so as to 
approximate ‘‘free-end’’ conditions. The initial out- 
of-roundness was measured and the tube was then 
placed in a creep chamber and the temperature was 
brought up to the desired level. The pressure was 
then gradually applied until it reached the precalcu- 
lated level and held steady until the tube collapsed. 
The time of collapse was noted and, in ten of the thirty 
tests, the increase of deformation with time was also 
measured. 


(4.2) Description of Test Apparatus 

A schematic drawing of the test setup is shown in 
Fig. 6. 

The temperature within the test chamber 1 was con- 
trolled and recorded by a Minneapolis-Honeywell 
recorder-controller 2. The temperature was con- 
trolled to within +3°F. at six points along the length 
of the specimen. Air was used as the pressurizing 
medium and was applied through the reservoir 3. The 
pressure was maintained constant (within about 1/2 
in. of water) by the regulator 4 at the air-inlet and the 
level of the water in the water column 5. Automatic 
control of the water level was achieved through a sys- 
tem of relays and solenoid valves as indicated in the 
figure. The time of collapse, evidenced by a sudden 
drop of the water column below a predetermined level, 
was also automatically recorded by an electric clock. 
In the tests where the deformation of the cylinder was 
continuously measured, the signal was fed into the 
oscillographic strain-recording equipment 6. 

Fig. 7 shows the details of the creep chamber. The 
chamber consists essentially of two concentric cylinders, 
the outer cylinder serving as an insulating enclosure, 
with the strip heaters mounted on the outside of the 
inner cylinder. The center post was made of steel 
and the two ends of the specimen were supported By 
steel discs, 6 in. in diameter and machined to within 
+().001 in., placed a fixed distance apart. The method 
of supporting the ends of the cylinder is shown sche- 
matically in Fig. 8. The distance between the end of 
the specimen and the disc was on the order of 1/8 in., 
to allow for the differential expansion between the 
steel post and the aluminum specimen. The ends were 
attached to the discs with two layers of flexible tape, 
one of Teflon and the other of silicon-glass cloth. This 
provided a near-perfect seal against pressure loss, and 
the tapes showed no deterioration up to the maximum 
temperature (500°F.) of the tests. 

Fig. 9 shows the center post with the high-tempera- 
ture LVDT’s used for measuring continuously the 
deformation of the cylinder. The LVDT’s were placed 
120° apart at each section and three sets of LVDT’s 
were used along the length of the cylinder. The figure 
also shows a mounted cylinder. Fig. 10 shows the 
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CREEP COLLAPSE OF 
facility for measuring the initial out-of-roundness of 
the test cylinder. The specimen was taped to the 
center post as shown, ready for the creep collapse test. 
The housing held the center post vertically, to +0.0002 
in., on ball bearings, while being rotated by a geared- 
down electric motor at approximately 1 rpm. A shaft 
supported and positioned the differential transformer, 
and the strain signal was fed into an oscillographic 
recorder so that a continuous plot of the signal around 
the cylinder was obtained. Zero calibration was taken 
at the 6.000-in. diameter (1) and the deviation calibra- 
tion was taken at the 6.040-in. diameter (2) on the eali- 
bration plate which was machined to an accuracy of 
+().0002 in. Phase orientation was maintained with 
a fine wire stretched the length of the cylinder, which 
caused a slight blip on the record. Measurements 
were taken at five sections along the cylinder marked 
A, B, C, D, and E, respectively. The method of re- 
ducing the recorder data to give a measure of the 
ellipticity to be introduced in the equations is outlined 
in Appendix A of reference 10. 


(4.3) Test Results 


Table 1 gives a summary of the cylinder data and 
test results. Table 2 gives values of 7 for aluminum 
alloy 3003 H-14, based on a limited number of creep 
tests on tensile specimens. 

The method of obtaining the mean radius shown in 
column 2, Table 1, is outlined in Appendix A of refer- 
ence 10, and is not reproduced here. The mean radius 
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TABLE 2 
Values of Material Constants for AL Alloy 3003 H-14 








°F E T 





300 9. 5x10° 89, 420 
325 9. 4x10° 23, 37 
350 9.27x10° 6, 646 
375 9. 14x10° 2,121 
400 9. Ox10° 672 
425 8.81x10° 236 
450 8. 6x10° 89.2 
475 8. 33x10° 35.2 
500 8.0x10° 14.9 
4H 

7 = 2.3026 s2g eRT 

T, constant, 714 psi 

E elastic tensile modulus, psi 

T temperature, °K 

s constant, 2. 3x10° per hour per °K 


SH ss activation energy, 35,000 calories per mo 


R gas constant, taken as 2 calorie; mol, *k 


Values of ¢9, s and AH tabulated above are approximate values, taker 


from limited creep data 





























TABLE 1 
Summary of Experimental Theoretical Results 
1 2 3 4 5 6 7 ~ 9 10° 
. tor te 

Cylinder Number a h Woo Temp q Gor Theor Actual 94cr/Gor% 
2+ 3.0368 0.0436 0.0021 500 3.95 6. 58 10.2 27.0 +18.1 
7 3.0288 0.0437 0.0034 325 7.42 7. 84 24.9 5.8 -0. 59 
9 3.0370 0.0438 0. 0060 500 2.89 6.67 17.5 46.5 +74.4 
10 3. 0352888 0.0438 0.0053 475 4.52 6.95 10.2 3.3 -5.7 
10A 3. 0298%88 0.0438 0.0242 475 3.16 7.01 12.8 22.8 +23.2 
12% 3.0354988 0.0413 0.0031 500 3.36 5. 60 9.0 29.0 +66. 1 
13 3.0310 0.0410 0.0028 450 4.14 5.90 27.0 5.7 -15.4 
13A 3.030488 0.0410 0.0624 425 3.33 $. 06 2.8 7 -3.7 
14 3.0350 0.0410 0.0044 400 5.21 6.16 25.9 22.8 -0.96 
14A 3. 0352008 0.0413 0.0808 500 2.2 5.57 1.4 2.5 +16.6 
1688 3. 0380 0.0411 0.0021 350 5.80 6. 33 57.4 20.8 -1.8 
16Ae* 3.0346 0.0411 0.0314 350 4.46 6.37 37.5 21.9 -2.2 
17 3.0443 0.0416 0.0044 500 4.25 5.67 2.3 1.0 -6.5 
18%8 3. 0323008 0.0413 0.0042 425 5.25 6.17 9.1 0.71 -4.8 
18A 3.0297 0.0413 0.0816 375 3. 20 6.4) 19} 4.4 -7.9 
19 3.0280 0.0420 0.0044 400 5. 68 6.67 217 1.5 -4.4 
19A*® 3.353088 0.0415 0.0720 450 2.45 6.12 a1 17 § +19.8 
23 3. 0340 0.0410 0.0017 450 3.61 5.90 6] 0 41.2 9.8 
23A 3.0247#88 0.0410 0.0943 475 1.99 5.76 4.6 ° 8 +29.0 
24 3.032698 0.0410 0.0021 450 5.02 5.90 7.3 3 -4.2 
25 3.0332ee0 0.0412 0.0015 425 5.22 6.14 23.7 0.75 -8.1 
26 3.0360 0.0410 0.0072 475 4.27 5.70 3.3 5.4 +6. 6 
27 3.0322 0.0409 0. 0064 350 5.87 6.32 3.5 12 -0.53 
27A 3.0300 0.0409 0.0358 350 4.75 6. 33 4.3 3.7 1.2 
288 3.0356 0.0433 0.0057 375 6. 60 7.37 10.3 4.7 -1.0 
308 3.0350 0.0411 0.0055 400 5.27 6.21 17.1 1.4 -4.0 
32 3.0300 0.0415 0.0086 375 5.85 6.52 28 8.2 +2.8 
3300 3.039888 0.0410 0.0027 425 5.10 6.00 15.2 0.2 6.6 
3408 3.0357%e8 0.0411 0.0030 475 4.88 5.74 2.2 0.2 6.0 
35 3.0327 0.0414 0.0053 425 4.99 6.23 14.8 14.4 0. 20 





A Denotes rerun cylinder. 
* Figures are inaccurate when 4,, is large 
** Denotes instrumented test. 


*** Denotes mean radius of 5 sections a= 1/6 (ay +apt zac +ap+ap). Others in column 2 are mean radii at center section C (see Fig. 10) 
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Fic. 11. Variation of modulus of elasticity EZ with temperature 
for aluminum alloy 3003 H-14. 


recorded in column 2 is generally the outer radius taken 
at the center section of the tube. Measurements 
showed that the radius varied along the length of the 
tube. The mean thickness / in column 3, Table 1, 
was obtained by measuring the thickness of the cylinder 
wall near the ends by means of an ordinary micrometer. 
No attempt was made to measure the thickness at the 
center section, and some variation is to be expected. 
The mean initial ellipticity was calculated as detailed 
in Appendix A of reference 10. The critical pressure in 
column 7 was obtained from the formula 


Yer = [E/4(1 — v?)](h/a)? (37) 


which is the theoretical elastic buckling pressure for an 
infinitely long and perfectly circular tube. The value 
of v? was taken as 0.10 while the value of E was ob- 
tained from Fig. 11, which is based on data furnished 
to the investigators by the Aluminum Company of 
America. 

In Table 1, the cylinders with the suffix A are reruns, 
that is, ‘‘collapsed’”’ cylinders that were tested again. 
When the ratio ¢/q-, was relatively small, it was found 
that the collapsed cylinder merely assumed a _ pro- 
nounced elliptical shape and appeared otherwise un- 
damaged. As one of the chief difficulties encountered 
in this investigation was the variation, both in magni- 
tude and orientation, of the maximum ellipticity along 
the length of the cylinder, it was felt that cylinders with 
a definite elliptical shape would provide an opportunity 
of reducing the uncertainty regarding ellipticity. As 
can be seen from the test results, a closer prediction 
of t., was possible in these cylinders. 

Column 10 in Table 1 shows the percentage error 
in g., that would have to be assumed in order for the 
theoretical ¢,, to be the same as the observed f,,, all 
other parameters being exactly known. The value of 
t., is, however, not nearly as sensitive to an error in 
foo as it is to an error in Ger. 

Fig. 14 shows a comparison of the experimental and 
theoretical results on a semilogarithmic plot. The 
authors are of the opinion that Fig. 14, while giving a 
visual picture of the degree of agreement between test 
results and theory, is not nearly as informative as a 
careful study of Table 1. 


Fig. 12 is a deflection-time plot of some of the instru- 
mented cylinders. The deflection is the radial deflec- 
tion of the cylinder wall. Fig. 13 shows the end 
view of some collapsed tubes. 


(5) Discussion and Conclusions 


It may be noted that while an order-of-magnitude 
agreement exists between test results and theoretical 
predictions, the accuracy of the theory leaves some- 
thing to be desired. As foreshadowed in the preceding 
paragraphs, it would perhaps be unduly optimistic to 
expect closer agreement on account of the numerous 
possible sources of error in the various parameters, and 
the simplifying assumptions that have been made. In 
order that future investigations may help in eliminating 
some of these errors, it is of importance to discuss 


them in detail. 


(5.1) Accuracy of q., 

The first question that may be raised is the accuracy 
of ger as estimated by Eq. (37). This formula applies 
to an infinitely long and perfectly circular tube. That 
the boundary conditions were indeed close to those of 
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“free ends” are amply borne out by the instrumented 
tests, which showed that the deflections at the ends and 
center section of the tube along the same axis were of 
the same order of magnitude for the duration of the 
test. Thus one may conclude that little accidental 
“fixity”? of the ends came into play. The assumption 
of infinite length therefore appears justified. The 
second assumption of Eq. (37), namely, that the tube 
was perfectly circular, is obviously not true. In order 
to estimate the effect of ellipticity the empirical formula 
of Southwell® (see p. 221 of reference 5), 


Oy.p. 


Ger = 


a Sy.p.(1 — v*) 4a° (38) 

a = 
h E h® 
where cy.,, is the yield point of the material, was used 
to estimate g.,.. Although in some instances this 
formula yielded a much closer estimate of the collapse 
time than did a prediction based on Eq. (37), it failed 
to do as well in other instances. It appeared therefore 
that some uncertainty regarding g., was unavoidable, 
and it was decided to use Eq. (37). 

Fig. 5 shows the error to be expected in the predicted 
value of ¢., for a small error in the calculated value of 
Jer = Eh*/(4(1 — v*)a*]. No special tests were run to 
determine the value of / in this investigation, and a 
5 per cent error in its value is not beyond reason. 
While a is known accurately enough at any given 
section, its value varied from section to section along the 
same tube. However, / is known only approximately, 
as the thickness was measured at several points along 
the circumference only at the ends of the cylinder, and 
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Fic. 13. Shape of cylinders at collapse. (Middle cylinder shows 
shape before test.) 


the arithmetic mean value wus taken as the true value 
of h. Here again a 5 per cent variation in / is well 
within the bounds of possibility. Thus a 15 per cent 
error in the value of h* is possible. In any event, it 
would appear that a 5 to 10 per cent error in the calcu- 
lated value of g., must be regarded as practically un- 
avoidable, even if very accurate data are available. 
As may be seen from Table 1, column 10, an error of 
this order is sufficient to account for much of the ob- 
served discrepancy between theory and tests. 

An attempt was made in the instrumented tests to 
obtain an estimate of q¢,, by using the Southwell plot 
(see p. 177, reference 5) of {/q against f as the pressure 
was increased to the test pressure. Such attempts 
turned out to be generally unfruitful, possibly for two 
reasons, (1) the change in magnitude and orientation 
of the ellipticity along the length of the cylinder and 
(2) the relatively high values of g/q., that are required 
for an accurate estimate of g., by the Southwell method. 
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Fic, 14. Comparison of tests with theory. 
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(5.2) Accuracy of fy 


Fig. 5 shows the error to be expected in the predicted 
value of ¢,, as a result of a small inaccuracy in fo. It is 
to be remembered that the inaccuracy in foo arises not 
merely from an error in measurement of wo) but also 
due to its variation over the length of the cylinder. 
Unless it is possible to test cylinders which have the 
same out-of-roundness all along the length, and along 
the same axis, this error is unavoidable. As foo is 
proportional to w/h, an error in h would also affect 
foo, however accurately wo is known. As a matter of 
fact woo is known, relatively, with sufficient accuracy 
(up to +0.0002 in.) in all the cylinders, at any given 


section. 


It was observed that whatever the original orientation 
of the ellipticity of the cylinder at the beginning of the 
test, the tube finally assumed an elliptical shape with- 
out change of direction of major and minor axes along 
its whole length. It is evident that in cylinders with a 
pronounced “twist,” this involves a rotation of the 
axes as well as a deflection of the cylinder walls during 
the creep process. Evidence of this was forthcoming 
in several of the instrumented tests. Some deflectom- 
eters actually showed a reversal in the direction of 
deflection with time. In the case of the reruns, this 
additional complication did not take place. The bet- 
ter agreement between theoretical and actual collapse 
times in the case of reruns supports this conclusion, 
although here again the magnitude of the ellipticity 
varied along the length of the cylinder, although not 
its orientation. 


A word of explanation seems appropriate with regard 
to Figs. 12(a) and (b) which show a plot of the radial 
creep deformation (that is, the total deformation minus 
the instantaneous elastic deformation) against time 
recorded in several of the instrumented cylinders in 
which the defiectometers for measuring the radial de- 
flection were mounted as shown in Fig. 9. In Fig. 
12(a), the deflectometer 4 at the center section of 
cylinder 28 shows a reversal in the direction of de- 
flection, a phenomenon which has been mentioned 
previously. This seems to indicate a relative rotation 
of the different sections of the cylinder before it assumed 
its final, more or less uniform, elliptical shape at collapse. 


The deflections at the top, center, and bottom sec- 
tions of cylinder 16 are also shown in the same figure. 
The deflections at these sections are of the same order 
of magnitude as would be expected in a cylinder of in- 
finite length. It may be concluded, therefore, that 
the tape seals at the ends imposed little or no accidental 
restraint on the end sections of the cylinder. The 
same conclusion may be drawn by an inspection of the 
curves for cylinder 30 shown in Fig. 12(b). 


The rapid increase in the deflection as the cylinder 
approached ‘‘collapse’”’ is clearly evident in all the 
curves and particularly in those (like cylinders 18 and 
33) which collapsed in a very short time. 
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(5.3) Inaccuracy in Estimating Collapse Time 
“Collapse” of a cylinder is a somewhat inaccurate 
word. The test setup was such that “‘collapse’’ was 
synonymous with a deflection of the cylinder wall 
of about 1/8 in., i.e., of the order of magnitude of three 
times its wall thickness. This appeared a reasonably 
good measure of collapse in the sense that the change 
in shape of the cylinder was such as might render it 
functionally useless. In any event, much larger de- 
flections are meaningless from the point of view of 
linear shell theory. However, it must be remarked 
that the instrumented tests showed conclusively that, 
near collapse, the rate of deformation increased very 
rapidly, and, presumably, the cylinder would have 
deformed considerably more if the tape seal at the 
ends of the cylinder did not give way. Nevertheless, 
some uncertainty must exist regarding the “‘time of 
collapse,’’ probably about 10 to 15 per cent of the 
measured /¢,,.. In other words, a nonlinear shell theory 
could conceivably yield a different critical time at 


large deformations. 


(5.4) Accuracy of Material Constants 

The constants in the creep law were derived from a 
very limited number of tests, and as such are subject 
to the usual errors due to experimental scatter and to 
the variation of the constants from lot to lot of the 
material. The specimens were in the as-drawn condi- 
tion and some anisotropy is to be expected although its 
effect was probably small. 


(5.5) Accuracy of Creep Law 

The creep law used in deriving Eq. (20) does not 
take into account the ‘“‘transient’’ stage of creep. In- 
clusion of this stage results in a considerably more 
complicated equation and, in view of the other in- 
accuracies, did not seem warranted at the present 
stage of the development of the theory. A somewhat 
simple semiempirical creep law, including the transient 
stage of creep, has been proposed by Pao and Marin‘ 
but, being essentially a curve-fitting procedure, it re- 
sults in several empirical approximations which make 


it unsatisfactory. 


(5.6) Assumption of a Sandwich Wall 

A quantitative measure of the inaccuracy resulting 
from this approximation is difficult but its magnitude 
is believed to be small. It should be possible to 
solve the equations assuming a solid wall by integrating 
through the wall thickness, an obvious problem for 
further research. 
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Some Aspects of Nonequilibrium Flows 


RAYMOND SEDNEY* 
Ballistic Research Laboratories 


Summary 


In this paper are discussed some of the general features of non- 
equilibrium flow. In particular, vibrational relaxation is dis 
cussed in detail. This case is somewhat simpler than dissociation 
and ionization but it illustrates some of the main new features of 
nonequilibrium flow. Those aspects of two-dimensional and 
axisymmetric flow behind shock waves are examined analytically 
which yield significant information without requiring numerical 
solution of the governing equations. 

The thermodynamics of a vibrational relaxing gas are dis- 
cussed. The conditions for simulating flows are noted. Crocco’s 
theorem and the characteristic equations are derived. Then a 
simple method of obtaining the initial gradients of the flow vari- 
ables behind a shock is shown. These gradients are used in dis- 
cussing two particular flows. An exact solution for flow over a 
cusped body is obtained. Flow over a wedge near the tip and 
far from the tip is considered. It is found that far from the tip 
a boundary-layer type phenomenon occurs 


Symbols 
“te pressure coefficient 
Geel specific heats at constant pressure and volume, 


respectively 


D 1 — M,? + cot?A 
Dy = operator for differentiation along characteristics 
E = internal energy 
é = [Ei(T.) — Ei(T.)|/RTa 
Fi 2.3.4 coefficients in gradient functions 
f = [(Ei(T.) — Ei(T,)|/¢paTarq 
fo = JT, 0S/0s 
h = enthalpy 
K = curvature of shock wave 
K = curvature of streamline 
l = displacement of shock, Fig. 7(c) 
VW = Mach number 
n = coordinate along normal to streamline 
p = pressure 
qd = velocity 
R = gasconstant 
r = radial polar coordinates 
S = entropy 
s = coordinate along streamline 
r = temperature 
t = time 
£ = coordinate in free-stream direction 
y = coordinate normal to free-stream direction 
e =0-¢ 
B = shock-wave angle 
= ratio of specific heats 
€ = 0 or 1 for two-dimensional or axisymmetric flow, 


respectively 


¢ = vorticity 

” = coordinate along shock wave 

6 = direction of velocity vector 

6 = characteristic vibrational temperature 


Presented at the Hypersonic Flow Theory Session, IAS 28th 
Annual Meeting, New York, January 25-27, 1960. 
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r = B — @, behind shock wave 
” = sin-!(1/M,) 
é coordinate normal to shock wave 
é, relaxation distance |Eq. (5.12)] 
p - density 
o = distance along shock 
T = relaxation time 
o = angular polar coordinate 
y 8-60 
w vorticity 
Subscripts 
a active degrees of freedom (translation and rotation) 
i = internal degree of freedom (vibration ) 
e = equilibrium conditions 
f = frozen conditions 
t = stagnation conditions 


free-stream conditions 


(1) Introduction 


8 tere HAVE BEEN MANY recent publications con- 
sidering the flow of a gas not in thermodynamic 
equilibrium. No complete bibliography will be given 
here, but several papers pertinent to the present work 
will be mentioned. 

Probably the first work on this subject, for tine case 
of shock waves in the flow, is that of Bethe and Teller.' 
A qualitative discussion of flow over a wedge was given 
by Ivey and Cline? in which the departure from thermo- 
dynamic equilibrium was due to vibrational relaxation 
of a diatomic gas. (Some quantitative results will be 
given here for this case.) An approximation to the 
hypersonic flow over a sphere, where the departure 
from equilibrium is due to dissociation, was considered 
by Freeman.* Flows with small disturbances, in par- 
ticular, the flow over a thin wedge, were considered by 
Moore and Gibson.* 

In a pure diatomic gas at reasonable temperatures, 
three processes can cause departure from equilib- 
rium: vibrational excitation, dissociation, and ioniza- 
tion. If the relaxation times for these processes are 
sufficiently different, they may be treated separately. 
The magnitudes of the effects on the flow variables due 
to the three processes are quite different. The gross 
effects are directly related to the energy necessary to 
excite the new degrees of freedom. 

In this paper vibrational relaxation only is considered. 
This is somewhat simpler to handle than dissociation 
and ionization, but should illustrate some of the main 
features of nonequilibrium flow. Because the smallest 
nonequilibrium effects result from vibrational excita- 
tion, some of the results obtained here differ little from 
the results of calculations based on equilibrium flow, 
or differences are significant only in a negligibly small 
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region of flow. It is hoped that the methods used here 
can be applied, with almost equal ease, to other non- 
equilibrium processes. 

Only certain aspects of two-dimensional and axi- 
symmetric flows are examined in this paper; these will 
not involve numerical solution of the differential equa- 
tions. After a discussion of the appropriate thermo- 
dynamics of a vibrational relaxing diatomic gas, the 
governing equations are transformed to yield the 
generalization of Crocco’s theorem which relates the 
entropy change to the vorticity. Next a derivation 
and discussion of the characteristic equations are given. 
Using the appropriate shock transition relations and 
the equations of motion, the streamwise gradients of 
the flow variables are obtained employing natural 
coordinates (as indicated by Sternberg® for equilibrium 
flows). These prove to be quite useful in discussing 
two particular flows; namely, two-dimensional flow 
over (i) a cusped body which supports a straight shock 

yave, and (ii) a wedge. 

For (i) an exact solution is obtained that requires 
only two quadratures. For (ii) the flow is examined 
near the tip and far from the tip. At the tip the shock 
curvature and the wedge pressure gradient can be 
obtained rather simply by use of the gradient func- 
tions. Far from the tip the flow must be divided into 
two regions. To a first approximation a large region 
of the flow is the equilibrium wedge flow, but there is a 
small region near the shock wave where relaxation is 
important. Mathematically this small region ex- 
hibits a boundary-layer type of phenomenon. An 
indication is given of a means of obtaining the next 


approximation. 


(2) Thermodynamics 


Basically, the difference in the descriptions of equilib- 
rium and nonequilibrium flows results from differences 
in the thermodynamic behavior of the gas, the dy- 
namic aspects being the same. A clear explanation of 
an appropriate model of a nonequilibrium system (and 
the necessary assumptions) was given by Wood and 
Kirkwood. The same model will be used here. 

For the case of a gas subject to vibrational relaxa- 
tion, the degrees of freedom are divided into two classes: 
the active (translation and rotation), for which the 
subscript @ will be used, and the internal or inert (vibra- 
tion), for which the subscript 7 will be used. It is 
assumed that local thermodynamic equilibrium exists 
within the classes but not between them. The rate at 
which equilibrium is approached is governed by a rate 
equation which can be specified in terms of the state 
variables of each class. There are other assumptions, 
to be mentioned later, which are convenient to make. 
The model presented® is intermediate between a macro- 
scopic and a microscopic description. 

From the assumption that the Gibbs relation holds 
for the “a” class, it follows that the entropy change 
of the ‘“‘a”’ class is 


dSq = Cp, dT ,/T, — Rdp/p 
where the perfect gas law is assumed, 


pb = pRT, (2.1) 


It is assumed that the ‘‘z”’ class is specified by its tem- 


perature, 7;, so that the entropy change is 
dS, = dE,/T; 
The total entropy change is then 
dS = dS, + dS, 

= Cp, d1,/T, — Rdp/p + dE;/T, (2.2) 

The rate equation assumed is the linear form 
dE,/dt = [E,(T.) — E,)/r (2.3) 

where 7 is the relaxation time and F£,(7,) is the value 
of /; if equilibrium existed at the temperature 7,,. 
For a flow the derivative in (2.3) is the substantial 
derivative. 

The linear form of the rate equation presupposes that 
the departure from equilibrium is not too great. It is 
exact if the vibrators are quantized harmonic oscilla- 
tors.' The functional relation F/,(7,), for various 
gases, can be obtained from tables of the properties of 
gases.’ With the assumption of harmonic oscillators, 
E,(T,) can be calculated from 

E,(T,) = R@,/{exp (0,/T.) — 1] (2.4) 
where 0, is the characteristic vibrational temperature. 
This form will be used when some specific flows are 
calculated. It is convenient, though not necessary, to 
carry along both variables E; and 7;. The functional 
dependence E;(7;), for the case of harmonic oscillators, 
has the same form as /(7,) in (2.4). 

For future reference the following well-known rela- 
tions are recorded. For equilibrium flow, 7, = 7; and 


Co —C,=R 
= Cpa + Ci 
c, = dE,/dT, 


where Cp 
Cy = Coq + Ciy 


Also, Cog ~— Cog = R 


(3) Equations of Motion 


In this section the equations of motion of a pure 
diatomic gas subject to vibrational relaxation will be 
discussed. The equation of state (2.1) is assumed. 
By a minor modification the results can be extended to 
the case of a mixture of diatomic gases where only one 
constituent relaxes. It is necessary to know the effect 
of the other constituents on this one. This case would 
arise, in practice, if the relaxation time of one constit- 
uent were much shorter than that of the others. This 
would be the case for air since the relaxation time of 
oxygen is about one-fifteenth that of nitrogen. This 
extension, however, will not be included here. 

The equations will be written in natural coordinates 
s,n, where s and n are distances along the streamlines 
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and their orthogonal trajectories, respectively. The 
form of the equations of conservation of momentum, 
mass, and energy do not change because of relaxation; 


they are 
pg Og/Os = —Op/ds (3.1) 
pq? 00/0s = —Op/On (3.2) 
1 Op 1 Og 06 e sin 6 at (3.3) 
p OS qg Os On y 
h + (q?/2) = h, (3.4) 
where g is velocity, @ is the velocity direction, « = 0 


or | for two-dimensional or axisymmetric flow, respec- 
tively; y is the coordinate perpendicular to the axis of 
symmetry; h = c,, 7, + E; + p/p is the enthalpy; and 
h, is the stagnation enthalpy (constant along stream- 
lines). These four equations, together with the equa- 
tion of state (2.1) and the rate equation, give a com- 
plete description of the flow. In natural coordinates 
the rate equation (2.3) is 


gq 0E,/0s = [E(T,) — E,]/r (3.5) 


The relaxation time 7 is a function of pressure and 
temperature. In the following it is assumed to be con- 
stant. This would be a poor assumption if there are 
large changes in pressure and temperature; however, in 
the regions of flow considered changes of appreciable 
magnitude in these quantities do not occur. 

For the above six equations, the conditions for similar 
flows over geometrically similar bodies are easily de- 
duced. In addition to the requirements of fixed Mach 
number and ratio of specific heats, rU/// must also be 
fixed for similar flows. Here U and / are representative 
velocity and length, respectively. This type of di- 
mensionless ratio was used by Freeman.* Experi- 
mentally, the appearance of the new length scale 7rU 
makes simulation more difficult. Analytically it means 
that the simple ‘‘similarity solutions’ (wedge flow and 
Prandtl-Meyer expansion) no longer exist. 

A useful relation for equilibrium flows is Crocco’s 
theorem.’ It relates the entropy gradient to the vor- 
ticity. For a relaxing gas this relation can be derived 
as follows. The expression for the entropy change 
(2.2) can also be written 


T, dS = dh — dp/p — [1 — (T,/T)] dE, 
Using (3.4) to eliminate h, 
T, dS = dh, — qdq — dp/p — [1 — (T./T))] dE, 


From this the directional derivatives of S in the s and 
n directions are evaluated. After elimination of the 
pressure by use of the momentum equations (3.1) and 
(3.2), these derivatives are expressible as 


T, 0S/ds = —[1 — (T,/T))] dE;,/ds (3.6) 


T, 0S/On = gt — [1 _ (T./T:)| OF ;/on + 
dh,/dn (3.7) 


where ¢ is the vorticity, 


¢ = q 00/0s — Og/On 


Note that 0h,/O0s = 0, but Oh,/O0n = O only for iso- 
energetic flows. If the flow is in equilibrium (7, = 
7), Eq. (3.6) gives the result of constant entropy along 
streamlines, and (3.7) reduces to the usual form of 
Crocco’s theorem. For three-dimensional unsteady 
flow the vector form of Crocco’s theorem is 


0q/0t —q Xw = T,VS + Vh, + [1 — (T./T)] VE, 


where w is the vorticity vector and V is the gradient 
operator. An equivalent result, but in different ther- 
modynamic variables, was given by Broer." 

It is convenient to eliminate p and p from the equa- 
tions of motion and introduce S. The resulting equa- 
tions for isoenergetic flow are 


. O”g 06 


(M2 — 1) _ me. ( q Fe (3.8) 
i Os q on y Cygl a] OS ai 


q? 00/0s — gq 0g/On = 
T, 0S/On + [1 — (T,/T))] OE,/On (3.9) 
T, 0S/ds + [1 — (7,/T)] OE;/Os = 0 (3.10) 
qg OE;/ds = [E(T.) — E|] ‘tr (3.11) 
Cogl'a + Ex + (qg?/2) = hy (3.12) 


where M,? = QCrq/CpghT 


Thus only the ‘‘frozen’’ Mach number enters. Note 
that (3.9) and (3.10) are simply (3.6) and (3.7), while 
(3.8) comes from (2.1), (3.1), (3.3), and the expression 
for 0S/0s. The unknowns are gq, 6, 7,, E;, and S. 

It is easily verified by the standard technique that 
the above partial differential equations are hyperbolic 
if 1, > 1. Solving the determinant equation for the 
four characteristic directions, one obtains the stream- 
line direction counted twice and the Mach lines based 
on M,; that is, the Mach angle is » = sin~! (1/M,). 
The equations of motion in characteristic form are 
(3.10), (3.11), and the following two equations: 


cot wDiq F gDs6 = 
sin pu [E(T2) _ E,] 
y lige 
cot u{ TDsS + [1 — (T./T)] DsE:} /¢™. (3.13) 


eg sin 6 sin wu 


where the operators D, are defined by 
sec uD, = 0/O0s + tan pO/On 


that is, D, indicates differentiation along the Mach 
lines inclined at + to the streamline. The appear- 
ance of the ‘‘frozen’’ Mach lines has been noted by 
other authors.* '® A numerical calculation, using 
the characteristic form, for flow over a wedge is in 
progress and will be reported in another paper. 


(4) Gradients Behind a Shock Wave 


A general method for computing the flow variable 
gradients behind a shock wave in a two-dimensional 
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SHOCK 






STREAMLINE 








— a 


Fic. 1. Natural coordinate system and notation for flow behind 
a shock wave. 


flow was developed by Thomas"! for a gas not subject 
to relaxation effects. A simpler method, using natural 
coordinates, was indicated by Sternberg.® For two- 
dimensional flow the gradients are particularly useful; 
e.g., the slope of the streamlines at shock polar points 
in the hodograph plane can be obtained from them (the 
so-called ‘‘hedge-hog”’ introduced by Busemann). It 
is found that the gradients are proportional to the shock 
wave curvature. 

For axisymmetric flow the same gradients can be 
computed, but they are linear combinations of the two 
curvatures: KA, in the meridional and 1/y in the 
azimuthal planes. This fact limits the usefulness of 
the gradients for this case. Extensive tabulations have 
been made recently’” of the coefficients of K, and 1/y; 
with these the gradients can be computed easily. (The 
two-dimensional result is obtained by setting the co- 
efficient of 1/y equal to zero.) 

For a uniform flow which is in equilibrium in front 
of a shock wave but not behind it, the gradients will 
depend on the relaxation time. Behind the shock the 
flow variables are computed from the shock relations 
for constant specific heats, i.e., from the frozen rela- 
tions.’ E; does not change across a shock.! To ob- 
tain expressions for the gradients, the equations of 
motion are combined with the shock relations; the 
normal derivatives are expressed in terms of deriva- 
tives along the shock wave and along the streamline. 
Referring to Fig. 1, 


0/Oc = cos \ 0/Os + sin \ 0/On 
where o is distance along the shock wave, 8 is the shock 
angle,and\ = B— 86. Since 
0/00 = 08/00 (0/08) = K, 0/08 
the normal derivative can be written 
0/dn = [K, 0/08 — cos 0/ds]/sin dX (4.1) 
The normal derivatives in (3.8) and (3.9) are eliminated 


by use of (4.1); 0£,/08 = 0, and 0q/08, 00/08, etc., 


can be computed from the shock relations. The fol- 
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lowing two equations then result: 


(M,? — 1) 0q/qds + cot \ 06/ds = 
(e/y) sin 0+ fi + (Kyk2/sin X) (4.2) 


cot \ 0g/0s + ¢ 00/ds = 
Kw (fe + gki)/(q sin \) (4.3) 


where k, = 0g/08 ko = 006/08 
f, = (OF,/O0s)/cp,T. = [Ei(Ts) — E,(T.)] /(cp,T a7) 
fo =e ip OS reyes 


and all quantities are evaluated immediately behind 
the shock wave. Solving (4.2) and (4.3) for 06/0s and 
0qg/0s, one obtains 


00/Os = FiKy + (€/y)F2 + (fi cot A)/D (4.4) 
(1/q.,) Og/Os = F3K,, + (€/y) Fy —gfi/(g.P) (4.5) 
where 
i = id — M,”)(Rki + fog!) + qk2 cot \] (gD sin d) 
Fy 
F; 
F, = —(q sin 6)/q..D 
D=1-M,?+ cot?A 


= (sin @ cot A)/D 


[(ki + fog") cot X — ghk2| /(g.D sin d) 


I 


The F; are functions of M/.., 8, and y, and have been 
tabulated’? for y = 1.4, 11 < M. < 10, and 
ain? (1/ MM.) < 8 S «/2. 


By use of (3.1), (3.3), and (3.4) the gradients of £, p, 
and 7, can be obtained. Since /; and D are positive, 
the effect of relaxation is to increase the streamline 
curvature 0@/Os and decrease the velocity gradient for 
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Fic. 2. Nondimensional streamline curvature behind a straight 
shock. 
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given Ky and y. For reasonably strong shocks at 
ordinary 7.., E,(7.) can be neglected compared to 
E(T.). 

Even for two-dimensional flows it is seen that the 
gradients will not be simply proportional to K,,; this 
will limit their usefulness. Some application will be 
made of these gradients in the next two chapters. 


(5) An Exact Solution 


As mentioned previously, the solution for flow 
over a wedge is no longer simple when the gas behind 
the shock is relaxing; the flow behind the shock is not 
uniform and the shock is not straight. This flow will 
be considered in the next section. Here it will be 
shown that an exact solution can be obtained for two- 
dimensional flow over a cusped body (with curvature 
becoming zero asymptotically) which supports a 
straight shock wave. 


From (4.4), with e = 0, it can be seen that if the 
shock is locally straight (K, = 0), the curvature of 
the streamline at the shock can easily be calculated. 
Letting K, = 00/0s, a convenient nondimensional 
measure of streamline curvature is the following 
quantity: 

K,7g~/e: = (Rg cot r)/cp,gD (5.1) 
where e, = [E(T,) — E,(T.)|/RT. (5.2) 


A plot of the expression in (5.1) is shown in Fig. 2 for 
a range of values of M., and 8. For the case of pure 
N. with 7, = 300°K, B = 60°, M.. = 6, where (2.4) 
was used with 6, = 3,336°K, it is found that 6 = 41.1 
and K, = 5.5 ft.-! The value of 7 was taken from 
reference 8. 

Since (5.1) has no physical length scale in it, the shock 
wave could be straight over any finite distance and K, 
would be the same everywhere. Considering the man- 
ner in which a numerical solution of the characteristic 
equations would proceed, starting from a straight shock, 
one can see that all the flow variables depend only on 
the distance along the normal to the shock. A “shock- 
oriented’”’ coordinate system £, is introduced, as 
shown in Fig. 3, and a solution with flow variables inde- 
pendent of » is sought. The streamlines will be parallel 
curves with initial curvature given by K, in (5.1) and 
initial slope calculated from the frozen shock relations. 
If the streamlines are physically meaningful, any one 
of them can be chosen as the body. Reflecting this 
curve through the free-stream direction gives a pointed, 
cusped body which supports a straight shock wave. 
Anticipating a result from the next section, the curva- 
ture of the body should approach zero and the angle of 
the body should approach the equilibrium wedge angle 
appropriate to the given shock angle. These expecta- 
tions are verified by the solution obtained below. 

Eqs. (3.8) through (3.12) are transformed to the 
(t,n) system by 


0/ds = sin (8 — 6) 0/d0E + cos (8 — 8) 0/On 
0/dn = —cos (8 — 6) 0/0E + sin (8 — 8) 0/On 


1) STREAMLINE 
BODY 





€ 


Fic. 3. Coordinate system for exact solution. 


Setting all n derivatives equal to zero and letting 
y = B-—6 


one obtains the following equations: 


: og 06 =gsinyOE 
(M,? — 1) sin + g cos some ' (5.3 
aed head al Eee 
q sin y 00/0 + cos ¥ 0g/d0E = O (5.4) 


T, 0S/dE + [1 — (T./T)] 0E,/doE = 0 (5.5) 
qsin y OF,/dt = [E,(T,) — E|] /t (5.6) 

Coala + Ei + (g?/2) = hy (5.7) 
On elimination of é, (5.4) can be written 

dq/q = —tan y dé 
Integrating, 

g cos Y = gq, cos (8 — Oy) (5.8) 
where the subscript f denotes values computed from 
the frozen shock relations. 

Combining (5.3) with the differentiated form of 
(5.7) and eliminating £, a first-order ordinary differ- 
ential equation is obtained. The result of integrating 
this is 
T. = [Ta, + (g7/R) sin? (8 — 6,)] X 

cot (8 — 6,) tan y — (g7/R) cos? (8 — 6,)tan?y (5.9) 


where again the frozen shock relations have been ap- 
plied. £; is obtained from (5.7), with the aid of (5.8) 
and (5.9): 
E, = kh, + [(Cpg R) tan? y — (1/2) sec? y] qr X 
cos® (8 — 6) — ¢p,[Ta, + (g7/R) sin? (8 — 6))] X 
cot (8 — 6, tan y (5.10) 


With all the variables expressed in terms of 6, (5.6) 
can now be used to obtain @ as a function of —&. Thus 


te] 
& = gsr cos (8B — a) F(6) dé 
6 


F(6) = tan y (dE,/d6)/[E(T,) — Ej] 
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Fic. 4. Shape of body which supports a straight shock wave. 





The equation of the body is, in parametric form (param- 
eter 0,), 


O% 
& = g,t cos (8 — a) f F(0) d6 
os 


Ob 
No» = Qyst Cos (B — a) f [F(6)/tan y] do 
Os 


The above integrations are carried out until E,(7,) 
is equal to F#;, that is, until equilibrium is attained. 
In principle, & and », become infinite for this condi- 
tion; actually, the curvature of the body approaches 
zero rather quickly and the gas is, for all practical 
purposes, in equilibrium at finite £, and np. 

A representative result (obtained by numerical 
integration) of the above formulas is shown in Fig. 4 
for No at M.. = 8, B = 65°, T.. = 300°K. Eg. (2.4) 
was used to obtain £,(7,), with 6, = 3,336°K. The 
angle of the body changes about three degrees from 6; = 
43.6° to 0, = 46.7°, the latter being the wedge angle 
that corresponds to a shock angle of 65° when the 
equilibrium shock relations are used. For the latter 
relations see, e.g., reference 13. The largest change in 
angle occurs when 6, and M., are close to the conditions 
for shock detachment in frozen flow. 

Thus an exact calculation (with two quadratures) of 
the flow with a straight shock wave is possible. These 
results can be used to check approximate methods of 
computing flows. In the next section this exact 
solution will be used to investigate a portion of the 
flow far from the tip of a wedge. 

It is desirable to define some measure of the distance 
over which relaxation effects are important. As in de- 
fining the thickness of a viscous boundary layer, there 
is some arbitrariness in any definition of a relaxation 
distance. Following Moore and Gibson,‘ one such 
definition would be 


f [1 — (E,/E,,,)] dé (5.11) 
0 


where £;,, is the final equilibrium value of F,. If 
E;.,. is replaced by E,(7,), the value of the integral 
would not be changed much, but the integration would 


be easier. Calling this relaxation distance £,, 


=| 61 — [E,/E(T.)] } dé 
0 


Oe 
= gyt cos (8 — 6,) f { | (dE;/d0) tan py] /E,(T,) | de 


(5.12) 


wer 


where 6, is the equilibrium wedge angle. For the par- 
ticular example cited above, 


£& = 0.136 g.7 = 0.300 gyr 


sr 


No reason was given for the particular definition of a 
relaxation distance in reference 4; it is analogous to the 
definition of displacement thickness of a viscous bound- 
ary layer. In the next section a different method of 
obtaining a distance is given, which results in a coeffi- 
cient of 7qg.. thirty per cent higher than the one given 
for &,. 


(6) Flow Over a Wedge 


Several years ago Ivey and Cline? discussed the two- 
dimensional flow over a wedge. The discussion was 
mainly qualitative; however, it was pointed out that 
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Fic. 5. Nondimensional shock-wave curvature at the tip of a 
wedge. 
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at the tip the shock wave should have the angle ap- 
propriate to the frozen shock relations, and far from the 
tip the (smaller) angle appropriate to the equilibrium 
shock relations.!* Also, the shock curve is concave 
downward. Here some further quantitative results for 
this flow are shown. 

It seems most natural to use polar coordinates (r,¢) 
with the origin at the tip and ¢ measured from the free- 
stream direction. Transforming Eqs. (3.8) through 
(3.12), with « = O, to polar coordinates gives 


(M,? — 1)(r cos a Og/Or + sin a 0g/0¢) + 
g(r sin a 00/0r — cos a 00/0¢) = 
(q/Cp,1a)(r cos a OE;/Or + sin a OF;/0¢) (6.1) 


@(r cos a 00/Or + sin a 00/0¢) + 
g(r sin a Og/Or — cos a 0g/0¢d) = 
—T,(r sin a 0S/dr — cos a 0S/0¢) — 
[1 — (1, T;)| (ry sin a OE;/Or — cos adE;/Od) (6.2) 


T, (r cos a OS/dr + sin a 0S/0g) = 
—|1-— (7, T)| (r cos a OE;/Or + sin a OF ;/0¢) 
(6.3) 


g(r cos a OF,/Or + sin a OF;/0¢) = 
(r/r)[E(Ta) — Ei] (6.4) 


Cole + E, + (q? 2) = h, (6.5) 


where a = 0 — @. 

To consider the flow near the tip a dimensionless 
variable, r/rq.., is introduced. Until the appropriate 
boundary conditions for wedge flow are introduced, 
the following discussion applies to any two-dimensional 
flow (e.g., the modifications to Prandtl-Meyer flow 
over a corner). Assume that all the flow variables 
can be expanded in a power series in r/rg.,; then it can 
be shown that, to within O(r/rq..) terms, the flow is 
frozen. The form of these series would be 


fir, &) = fold) + (r/rq) filo) +... 
The zero-order terms satisfy 
(M,,2 — 1) sin ap Ogo/O0G — Go COS ay OH/OG = 0 (6.6) 
COS ay Oqo/Od — go SiN ay O%/OP = O (6.7) 
OSo/0d = 0 (6.8) 
0E;,/0¢ = 0 (6.9) 
Eqs. (6.6) and (6.7) present two possibilities: either 
M,,? sin? a = 1 
which is appropriate for corner flow, or 
0g /Od = 06%/0d = O 


which, after applying the wedge boundary conditions, 
gives frozen flow over a wedge. To within an error 
O(r/rqg..), the shock wave is straight. 

The next step would be to solve the differential equa- 
tions for q, 6, etc. Although this is feasible, it will not 
be done because the most interesting information can 
be more easily obtained from the gradient functions 
of Section 4. Since, for the wedge, 00/0s = 0, the 
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Fic. 6. Pressure coefficient gradient at the tip of a wedge. 


curvature of the shock at the tip is, from (4.4), 
| = —(fi cot d) ‘DF, 


The frozen values of the variables are used in the right- 
hand side. In nondimensional form, 


Kytq./e; = —(Rq. cot d)/¢p,gDF, (6.10) 


A. plot of the expression in (6.10) is shown in Fig. 5 


for a range of M., and 8. For the case of pure Ne 
with 7.. = 300°K, 8 = 60°, M.. = 6, p. = latm., 


@ = 41.1°, K, = 4.1 ft.-! 


where (2.4) was used with 6, = 3,336°K; the value 
of 7 was taken from reference 8, and F; from reference 
12. For the same conditions with M. = 10, A,, is 
220 ft.—!, which means that the curvature would not 
be detectable in any ordinary scale of experiment. 
With K,, known, the initial velocity gradient along 
the wedge can be computed from (4.5), and then the 
pressure gradient from (3.1). A convenient non- 
dimensional form for the pressure coefficient is 


(OC,/Os)(7TGa/ei) = 
2(pq/ pode) ((R/CpD) — (FsKwrqo/e:)|} (6.11) 


where Cp = 2(— — Po)/ Poa” 


A plot of the expression in (6.11) is shown in Fig. 6. 
For the above example (M/.. = 6), OC,/Os = 16.0 ft.-, 

Thus the shock curvature and the rate of change 
of the flow variables on the wedge at its tip can be 
obtained rather easily; if the details of the flow are 
desired, the series expansion method should be used. 
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Fic. 7. Description of flow far from the tip of a wedge: (a) the 
first approximation giving the equilibrium region; (b) detail of 
transition region; (c) the two regions joined. 


Unfortunately, the same method, i.e., employing the 
gradient functions, cannot be used to investigate cone 
flow. The use of natural coordinates at the tip of a 
cone is impossible because of the nature of the singu- 
larity there. The series expansion method could be 
used, but the differential equations for the first-order 
terms would be much more difficult to solve than those 
of the wedge case. 

To consider the flow far from the tip, a series ex- 
pansion procedure is again adopted, but now 7q../r is 
the variable. The assumption that the flow variables 
are analytic in this variable will be seen to be not uni- 
formly valid. The difference between this case and the 
previous one can be seen by examining (6.4). When 
r/Tqo approaches zero, the right-hand side vanishes 
and nothing unusual happens. When 7q../r approaches 
zero, the left-hand side vanishes and the result is that 
EXT.) = E,; that is, the flow is in equilibrium. Since 
derivatives have been lost, a boundary-layer type of 
phenomenon (singular perturbation) must be expected. 
The order of the system (6.1)—(6.4) is reduced from 
four to three by this limiting process. Also the result 
that the flow is in equilibrium is inconsistent with the 
boundary condition £; = E;, at the shock. 

If the flow variables are expanded in a series of the 


form 
f(r.) = f(b) + (7G0/r) f(b) +... 


the zero-order terms for g and @ will satisfy equations 
of the same form as (6.6) and (6.7) except with zero 
superscripts instead of subscripts, and M, replaced by 
the equilibrium Mach number. Thus, to within 
an error 0(7q../r), the flow is equilibrium wedge flow. 
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The shock is straight and the flow is uniform and 
parallel to the wedge. 

It is helpful to keep in mind here the procedure for 
obtaining the viscous flow over a body: (i) the non- 
viscous flow is computed ignoring one boundary condi- 
tion (i.e., tangential velocity equals zero); (ii) the 
structure of the apparent discontinuity (vortex sheet) 
at the body is investigated using the boundary-layer 
approximation; (iii) the nonviscous flow over the body 
plus displacement thickness is computed, ete. For 
the present problem steps analogous to (i) and (ii) 
will be taken. Note that here part of the problem is to 
determine the boundary, that is, the shock wave. 

In Fig. 7(a) a streamline and the fictitious ‘“‘equilib- 
rium shock,” predicted by the first term of the series 
expansion, are shown in dashed lines. The transition 
across the shock is governed by the equilibrium shock 
relations. The angle of the equilibrium shock is fixed 
by M.. and the wedge angle. Extended, it must go 
through the tip in order to satisfy the conservation of 
mass. The actual shock (solid line) must be displaced 
upstream of the equilibrium shock, but its curvature 
must vanish to within an error O(7g./r). To this 
approximation the nonequilibrium flow behind a 
straight shock describes the flow near the shock but 
far from the tip [Fig. 7(b)]. This flow is obtained 
from the exact solution of Section 5. Thus the “‘bound- 
ary layer,” or rapid transition region, is determinable. 
Compared to the equilibrium region between the shock 
and the body this transition region is negligibly small. 
It and the shock merge into the equilibrium shock. 

The question remains how to join the two regions— 
or rather, how far is the shock displaced from the 
equilibrium shock? This is answered by again in- 
voking the conservation of mass. A transition region 
streamline must asymptotically approach the equilib- 
rium streamline that originates at the same point of 
the free-stream flow [see Fig. 7(b)]. The displacement 
of the shock from the equilibrium shock is denoted by 
/. From the condition defining / it is easily shown that 


1 = b tan B,/[tan B, — tan (8, — 6,)] (6.12) 


where £6, is the equilibrium shock angle, 6, is the wedge 
angle, and 6 is the distance (measured perpendicular 
to the shock) between the shock and the asymptote to 
the transition streamline [see Fig. 7(b) or 3]. There 
is no convenient expression for b, but it is easily deter- 
mined graphically when Eqs. (5.11) are plotted. For 
the conditions of the example of an exact solution 
given in Section 5, viz., 8. = 65°, 0, = 46.7°, M.. = 8, 
in No, (6.12) gives 


1 = 0.1777q.. = 0.0008 ft. 


whereas £, = 0.1369.7. In Fig. 7(c) the two regions 
are shown according to this approximation, i.e., to 
within an error 0(7qg./r). The equilibrium shock is 
shown extended to the tip. The next approximation, 
with an error 0(7q../r)?, would show the downstream 
edge of the transition region to be curved. 


(Continued on page 208) 
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Supersonic Panel Flutter of a Cylindrical Shell 
of Finite Length’ 


MAURICE HOLT* ano SAINSBURY L. STRACK** 


University of California at Berkeley and Boeing Airplane Company 


Summary 


flutter characteristics of cylindrical shells of finite 
In pre- 


Panel 
length in a uniform supersonic stream are determined. 
vious work, the aerodynamic term is calculated for a shell of 
infinite length and the equation of shell motion then reduces to 
a linear differential equation. In the present paper the aero- 
dynamic term is found for a finite cylinder and the corresponding 
shell equation is an integro-differential equation of fourth or 
eighth order. This is very difficult to solve as it stands but can 
be reduced to a differential equation by means of a plausible 
This 
simplified equation is solved here and results are compared with 
In general, it is found that the 


assumption about the behavior of the aerodynamic term. 


those of earlier calculations. 
use of a more realistic aerodynamic term leads to higher esti- 
mates of critical flutter Mach numbers. 


Symbols 
x, 7, 0 = cylindrical polar coordinates 
n number of circumferential modes of vibration 
lL Rh length, radius, and thickness of cylinder 
Ww radial distortion of cylinder 
U, M = free-stream velocity, Mach number 
M = critical Mach number 
B = (M? — 1)!/2 
¢ = disturbance velocity potential 
( = sound speed 
i = time 
x’ = x/l 
} = r/R 
i = Ut/R 
w = frequency of unsteady disturbance 
v = Rw/U, reduced frequency 
f, = reduced disturbance potential 
I Es = modified Bessel functions 
Un, Va, Wn = integrals of K,, 
Qa) = zeros of K,,'(z) 
/ = aerodynamic load on cylinder 
N,, No = applied axial and circumferential stresses 
d = (h/R)[12(1 — o%)]-/? 
D = bending stiffness 
E = Young’s modulus 
p, p = density of air and shell, respectively 


(1) Introduction 


I THIS PAPER we determine the conditions under 
which a cylindrical shell of finite length may de- 
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velop panel flutter in supersonic flow. The shell is 
placed in a uniform supersonic stream, so that, in its 
undistorted condition, its generators are parallel to the 
undisturbed stream direction. The shell is then as- 
sumed to distort in a periodic manner and small un- 
steady disturbances are introduced into the air stream. 
These disturbances change the aerodynamic forces, 
which in turn affect the mode of distortion of shell. 
When the aerodynamic disturbances reinforce the shell 
oscillations sufficiently to cause instability the shell 
is in panel flutter. We shall determine the critical 
conditions for panel flutter for different dimensions 
and modes of support of the shell and various undis- 
turbed free-stream Mach numbers. 

In previous treatments of panel flutter of cylindri- 
cal shells, notably by Miles,' Stepanov,” and Leonard 
and Hedgepeth,*® the aerodynamic forces on the dis- 
torted cylinder are calculated by a traveling wave 
approach. Such an approach is strictly applicable only 
to cylinders of infinite length or to one bay of an in- 
cylinder equally spaced ring 
stiffeners. This may partially account for the fact 
that flutter speeds predicted by the infinite span analy- 
sis are much lower than those observed in flight. In 


finite supported by 


the simpler problem of plane panel flutter, treated, 
for example, by Hedgepeth,* Fung,®° Luke and St. 
John,® and Eisley,’ finite dimensions, at least in the 
chordwise direction, are taken into account in aero- 
dynamic calculations and agreement between theoreti- 
cal and experimental results is fairly good. There 
would therefore seem to be a strong case for extending 
work on cylindrical shells to include the effect of finite 
cylinder length on aerodynamic forces. 

Leonard and Hedgepeth* and Miles! differ in deter- 
mining which circumferential mode of distortion is the 
least favorable for flutter conditions. In terms of the 
notation shown in Fig. 1, the mth such mode is defined 
by the term cos n#. Leonard and Hedgepeth’ claim 
that the cases n = 0 (uniform dilatation) and ” = 1 
(pure translation) do not represent panel motion and 


* 


Fic. 1. Geometry of circular cylinder. 
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Fic. 2. Variation of Mach number with thickness ratio re- 


quired to prevent flutter for an empty isotropic cylindrical shell. 
The results of Stepanov and Miles are for steel or aluminum. 
Those of Leonard and Hedgepeth are for aluminum. The stable 
region is to the left of each curve and the unstable region is to the 
right. 


consider only cases in which m > 2. The aerodynamic 
term is expressed in terms of cylindrical waves, used 
with Donnell’s equation’ for large values of 1 and with 
Fliigge’s equation’ for values of close to 2, when 
Donnell’s equation is believed to be less accurate 
(Hoff). Results are given in Fig. 2, showing the 
variation of critical Mach number J/, with ratio of 
shell thickness to shell radius. In the same paper 
Leonard and Hedgepeth® consider an infinite cylinder 
with regularly spaced ring stiffeners, but do not give 
any results in this case. 

Miles! specifically excludes the cases considered by 
Leonard and Hedgepeth (namely, » > 2) and pre- 
sents arguments to show that the fundamental insta- 
bility occurs at m = O for relatively short wavelengths 
(in terms of the length of the cylinder). To calculate 
the aerodynamic forces he assumes that the cylinder 
may be regarded locally as two dimensional and uses 
a linearized plane wave analysis to determine the load- 
ing. Donnell’s equation is used for the elastic motion 
of the cylinder and in this case it can be reduced from 
the usual equation of eighth order to one of fourth order. 
Miles also determines the negative damping ratio, which 
is a measure of the severity of the instability at any 
point. It is difficult to assess the significance of this 
quantity for a finite cylinder, since the traveling wave 
solution used in the infinite case is excluded by the end 
conditions in the finite case. 

Leonard and Hedgepeth and Miles have considered 
the effect of structural damping on the flutter bound- 
aries, and for the infinite shell they show that even a 
very small amount of damping may make the cylinder 
more prone to flutter. The explanation of this result 
is that the damping causes a phase shift in the wave 
motion, so that the air flow is able to feed more energy 
into the structure. For certain values of the param- 
eters this may cause a net gain of energy by the shell. 


In a recent paper, Miles'! investigates the effect of 
the boundary layer on the flutter boundaries and the 
negative damping ratio. Although no actual flutter 
boundaries are calculated Miles shows that the nega- 
tive damping ratio may be reduced by an order of 
magnitude, so that the flutter instability is less serious. 

A different approach to the infinite cylinder prob- 
lem is due to Stepanov.” In the first part of this paper, 
the aerodynamic forces are determined by piston 
theory!” and used both in Donnell’s equation’ and in a 
simple equation governing the shell motion due to 


Goldenveizer.'* Stepanov considers the cases n = () 
and m > 2, and finds that the critical speed has a mini- 
mum, for a fixed value of h/R (h = thickness, R 


radius of cylinder) when z = (), and has the same mini- 
mum for large values of 1 (of order 20). He obtains 
a value of critical flutter Mach number differing by 
unity from the corresponding value given by Miles 
(see Fig. 2). Stepanov uses a more direct approach 
than either Miles or Leonard and Hedgepeth and de- 
rives the stability boundary by solving an octic equa- 
tion instead of using a Nyquist diagram. Stepanov’s 
analysis could easily be extended to include structural 
damping effects and these would apparently not have 
the critical effect on flutter boundaries found in the 
other two methods of treatment. 

In the second part of his paper? Stepanov considers 
cylinders of finite length, subject to various end con- 
straints—for example, when simply supported at both 
ends, clamped at one end and free at the other, and so 
on. He again uses the simplified expression for the 
aerodynamic forces given by piston theory. To deter- 
mine the motion of the cylindrical shell he uses only 
the fourth-order equation given by Goldenveizer. 
The equation for eigenvalues is derived from a fourth- 
order determinant in this case, while Donnell’s equa- 
tion leads to a corresponding determinant of the 
eighth order. For m = 4 Stepanov obtains results, 
shown in Fig. 3, which, for most vaiues of the relevant 
parameters, give sufficiently high values of Mach 
number to justify the use of hypersonic small dis- 
turbance theory. However, the critical Mach number 
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depends on n~ in Stepanov’s analysis, and there is no 
reason why ” = 4 should be chosen in preference to 
n = 6, say, which would give much lower, nonhyper- 
sonic, values for /,. 

In the present paper, this analysis is improved by 
the introduction of more accurate expressions for the 
aerodynamic forces into the equations of shell motion 
in the hope that the results will be valid in supersonic 
as well as in hypersonic flow. The revised aerody- 
namic forces are used in conjunction with the Golden- 


veizer equation. 


(2) Formulation of Equations 


In this section we shall determine the aerodynamic 
forces on a distorted cylindrical shell, and also discuss 
the relative merits of different approximate equations 
governing the elastic response of thin cylinders. The 
aerodynamic term in these equations is the force on a 
circular cylinder of finite length slightly distorted in a 
general manner by unsteady effects; the derivation 
of this essentially follows Randall'* and Holt.’° The 
aerodynamic term is first determined from the linear- 
ized equation of unsteady cylindrical waves and then 
simplified by the static approximation in which the 
forces are calculated on a panel frozen in some mode at 
a particular instant. This approximation is justified 
in low frequency oscillations. The results may be 
applied to the cylindrical rear portion of a blunt-nosed 
hypersonic body provided that the nonuniform effects 
of the mixed flow region near the nose are ignored and 
the cylinder is treated as a duct, with an internal and 
external flow, which are independent of one another. 
The appropriate free-stream conditions in this case 
are the local conditions attained in the neighborhood 
of the cylindrical afterbody and not the conditions at 
infinity ahead of the blunt-nosed portion. 

Consider the cylindrical duct shown in Fig. 1. The 
x axis is in the direction of the undisturbed stream and 
cylindrical polar coordinates x, r, 6 are based on this 
axis. R is the radius of the cylinder in the undistorted 
state, w is the radial distortion, and w/R < 1. / is 
the length of the cylinder, and we shall in general take 
1/R ~ 4. The undisturbed stream moves with uni- 
form velocity U, Mach number M, and B = VM? — 1. 
We introduce a velocity potential 6 = Ux + ¢, where 
¢ is the disturbance potential which satisfies the equa- 
tion 
(1/r*) goa + ¢rr + (1 r)¢r pee B’ gz: “ 

2(M Cc) prt + (1 C7) put (2.1) 
and ¢ is the time. The equation of the distorted 
cylinder may be written 


r= R-+ w(x/I, 8, t) 22 


The following boundary conditions must be satisfied: 
In any meridian plane there is no disturbance ahead of 
the leading Mach lines through the rim of the duct, 
and the normal component of the velocity must vanish 
everywhere on the surface of the cylinder. The first 


condition implies that 
g=Q@ << 0 


g—~>Oasr— o~ forx> 0 


while the second gives 


Oy or 1 Ow 1 Ow 
= - + — (2.4) 
Or/,=R Ox l ax/I U oft 


We now introduce nondimensional coordinates 
x’ = (x/I)r’ = (r/R) = (Ut/R) 
Then the potential equation (2.1) becomes 
(1/r"*) G60 + grrr + (1/r') er — (RB/))*? 92 = 
2M?(R/D ori + Mg (2.5) 


and the boundary condition (2.4) 


(Ss) R ow 4 Ow (2 6) 
= — 2.0 
Or’ gaol l Ox’ ol 


We now assume that the potential varies periodically 
with time so that it is factored by the term e’” where 
v = Rw/U is the reduced frequency and w is the di- 
mensional frequency. 

Introducing the Laplace transform of ¢, 


0 


and dropping the primes, we obtain 


(1 r*) Dee + Orr + (1 1) Or =e (RB 1)? pe 
(R/1)-21vM*pe — M*v*?o (2.7) 


Since ¢ (and hence ¢) vary periodically with 6 we may 
write ¢ = f,(x, r) cos n#, and ¢ = f,(p, r) cos n6 where 


f, satisfies the equation 


Fine + C/r) he — 19? — (n®/r*)f fr = 0 (2.8) 
and) sq? = (RB/I1)? p? + 2ivM*p(R/1) — v?M? 


The boundary conditions for f, are 


f.(r) >0O as r>o 


. ole (2.9) 
(Of,,/Or),;=1 = [(Rp/l) + iv]@,$ 
where @(p, 0, t) = @,(p) cos nbe™ 


and the bar denotes the Laplace transform of the corre- 
sponding function in x. 
The general solution of Eq. (2.8) is 
f(r) = A,(p)K, (qr) + Br(p), (qr) 
(2.9) we take B, = 0, since 
f(r) = A,(p)K,(qr), and, 
to satisfy the second boundary condition (2.9), we have 


To satisfy the first of Eqs. ( 
I,—2> © asr—~ oo. Then 


q:An(p)Kn'(q) = [(Rp/)) + t)]@, 
_ [(Rp/)) + w]e, 


so that A,( 4 
so tha iAP) gK,,'(q) 
Rp, Ku(qr) 
and = {| — Dn° (2.10) 
inc f,(r) ( ] + i) qK,'(q) 
K,(qr) (2 =) 
If - ee 
VW n(Q, r) gK,,'(q) l r DP: 
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Fic 4. Contour described in Section (2). 
the inverse of Eq. (2.10) may be written 
x 
f,(x,r) = i) w,'(s)W,(x — 2, r)dz (2.11) 
0 


where JV,,(x, 7). the inverse of (a, r), is given by 
l i Balas i) 1 K,,(qr) 
W,, (x, 7) = | (| + ") : E e”* dx 
2m1 J c—io l b/ qk,'(q) 


The full solution for the mth circumferential mode is 


(2.12) 


then 


x 
o(x, 7, 0, t) = { w'(z)W,(x — 2, r)dz-cos nb-e'” 
0 
(2.13) 


In principle we can evaluate W,, if we replace p by q 
as the variable of integration, and use contour integra- 
tion. In practice this is complicated, and to simplify 
the problem we apply the static approximation. In 
Eq. (2.5) derivatives with respect to time are then 
omitted. This approximation is used by Hedgepeth'® 
in analyzing the flutter characteristics of flat panels of 
finite length. Hedgepeth refers to an earlier paper by 
Hayes” to find support for the approximation and, 
on this basis, states that it should be valid for Mach 
numbers greater than the value 1.6. Actually further 
justification of the approximation is desirable. This 
could be obtained by extending first Hedgepeth’s'® and 
then the present results. to take full account of time- 
dependent effects and comparing with existing results. 
The justification for using the approximation here is 
that it is a first-order approximation in terms of re- 
duced frequency v. Its introduction makes the analy- 
sis to follow fairly tractable, and gives a first estimate 
of the effect of introducing finite-length aerodynamics 
into the finite-shell flutter problem. 

Since we are only concerned with the aerodynamic 
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forces on the surface of the cylinder, and since w < R, 
we may evaluate the pressure on the undisturbed sur- 
face of the cylinder. Then on the cylinder surface, 


As, 1,0, t) = { w'(z)U,(x — 2, 1)dz cos nOe™ (2.14) 
0 
where 


a Jy wheels .. 

Unlx) = Une, 1) = 5 J ; x 

271 g-i@ 
K,(RBp/1) 


: e"dp (2.15) 
(RBp/1)-K,’(RBp/1) 


If f(x) is the inverse of f(p), (1/a)f(x/a) is the inverse 
of f(ap), so that the integrals to be evaluated are 


, l g+1 K,,( ) 
V, = f = P e”* dp 
251 J r~i pK,’ (p) 


where 1 is a positive integer and g is chosen so that the 
contour lies to the right of all the poles of the integrand. 

K,(p)/K,'(p) has a branch point at p = 0. The 
poles of the integrand occur only at the zeros of K,,’(p), 
and these zeros are known. They all lie to the left of 
the imaginary axis and occur in conjugate pairs. The 
number of zeros is equal to the even integer nearest 
to n + (1/2), giving no zeros for Ko’(p), two each for 
K,'(p) and K,'(p), four each for K3'(p) and K,4'(p), 
and so on. The required contour is shown in Fig. 4. 
lf K,'(p) has 2m zeros at p = —;(4 = 1,..., 2m), then 
the residues of the integrand are 

has nal | 
aiK,"(a;) 

For large p, A,(p) approaches zero sufficiently rapidly 
to ensure that the integral around a large semicircle 
I approaches zero as |p| ~ ~. Similarly, as |p| > 
0 the integral around a small circle y approaches zero. 
After some reduction it is shown by Randall'* that 


2m »aix 


»>—Pr 
ai ¢ 
y, = 2 +(-p" f x 
i=1 (a;? + n?) Jo p 


l 
214 
[K,’(p) 2? + 2 {L,"(p) 7 oe 
The functions V(x) are tabulated by Randall,'* and 
for the values n = 0, n = 2,” = 6 are shown in Fig. 
5. The integral term in Eq. (2.16) is fairly small, 
amounting to approximately 10 per cent of the total 
for nm = 2 and to a smaller percentage for larger values 
of n. 

The disturbance pressure, or aerodynamic load, on 

the cylinder is given by 


pU? de pU? , a i 
= . = - cos n0-e”” w'(z 
li“za as’ wna * 


. l 
Vi( (xe — 2) = Iz (2.17) 
(4 ) aa) (2.09 


where V,,(x) is defined by Eq. (2.16). Eg. (2.18) can 
also be written 
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‘ pt y " , ts 
Z cos nhe™ | w'(x) V,(0) + 


Bl 
= (@)- Va ( .) ds] @.8) 
WwW ° n ‘2. “= Zp Z eit 

BR Jo BR) ‘ 


We now formulate the elastic equations of the cylin- 


2 


der in a suitable approximate form. Many equations 
have been proposed, for example, by Moe'* and by 
Morley.'* The latter appears to be as simple as 
Donnell‘s equation, yet applies over a wider range 
of values of m. Kempner” and Hoff!® have com- 
pared Donnell’s equation’ with Fliigge’s equation,’ 
which is generally considered to be very accurate. 
These comparisons show that Donnell’s equation 
gives a very good approximation for large values 
of n but is less satisfactory for values of mn close 
to 3. For nm = 0 it again gives a very good repre- 
sentation of the displacement of a loaded cylindrical 
tube. Donnell’s equation of motion for a thin cylin- 
drical shell may be written in the form 


1 1 Ofw h 


Sq ° . — V3 
* x £ oF D ” 
ow  N, ew ) Z 
iV ° ‘ — Ps = V* 
Ox? R® 06° of? D 


where the coordinates are in dimensional form. The 
cylinder is of length /, radius R, and thickness h; N, 
and N, represent applied axial and circumferential 


stresses (positive in tension) ; 
a h\? l Eh* 
d? = ~ and D= a 
R/ 12(1 — o?) 12(1 — o?) 


the bending stiffness, where / is Young’s modulus and 
a Poisson’s ratio for the shell material; p, is the density 
of the shell and Z, as before, is the applied radial (pres- 
sure) load. We shall only use this equation in the 
case n = 0), taking NV, = N, = O, in general. In this 
case, in terms of the nondimensional coordinates pre- 


viously defined we have 


Ow ft UY? O*w l\4 1 i4 
-+ ph + R -—-w= Z (2.19) 


nut  D RR” oar 


0.68 
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Fic. 5. Examples of the V, functions calculated by Randall. 


006; 


002 — Total strain energy 


Bending energy : ae Stretching energy 
Se ae sat in ah a e_ 
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Number of circumferential modes, n 


Fic. 6. Strain energy due to bending and stretching. A/R 
0.01,1/R = O.8. 


For the cases nm > 2 we shall use Goldenveizer’'s 
equation'*? which is also used by Stepanov. Apart 
from the original derivation, we may justify the use of 
this equation by a study of Stepanov’s results. Stepa- 
nov first solves the panel flutter problem for a finite 
cylindrical shell by using the exact shell equation, and 
then also by a Galerkin procedure, which is found to 
converge well when only three terms are retained. 
Stepanov concludes from this result that the Galerkin 
approach gives reasonable convergence—a result also 
verified by Fung®—and applies it to the full eighth- 
order Donnell equations. The results converge towards 
the exact results obtained for Goldenveizer’s equation, 
and Stepanov states that for practical purposes, the 
use of the simple equations (i.e., Goldenveizer’s equa- 
tion) of shells of medium length is entirely justified 
when investigating problems of self-induced vibration 
of cylindrical shells. 


This equation is written, in the same notation as 
above, as 


Ow l\4 Oo” 2 O'w l\4 
, ad? l 
i (;,) ' (= r ) oot * (;.) x 


R°p, O'w (;.) R? O'Z 
(2.20) 


E 02 \R/ Eh 20! 


\< 


Before solving these equations it is convenient to 
discuss one further aspect of the problem—namely, 
the configuration of the cylinder in vibration. The 
possible modes of vibration of a cylindrical shell are of 
course infinite, and it would be very helpful if we could 
limit consideration to a finite range. The axial modes 
are determined by the boundary (end) conditions, but 
the circumferential modes (those corresponding to 
terms in cos ”@) are not. At first sight, therefore, un- 
less longitudinal stringers are applied as constraints 
to fix the value of m there appears to be no physical 
reason why the case » = 2 should be more important 
than the case m = 1,002. However, Arnold and War- 
burton?! give some insight into the way in which 
might be restricted by strain energy considerations. 
They consider the free vibrations of a simply sup- 
ported, thin cylinder of finite length, and determine 
the strain energy involved in vibration, for various 
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values of m. Further, for each value of n they also 
calculate the proportion of strain energy due to bend- 
ing and that due to stretching. Their results are shown 
in Fig. 6, where it is seen that as m increases more 
energy goes into bending, and less into stretching. 
Their results clearly show that there is a value of n 
equal to 6 or 7 (m must be an integer) for which the 
strain energy is a minimum. It is plausible to assume 
that the shell will vibrate most readily in the con- 
figuration requiring minimum strain energy. If this 
assumption is correct, then we may take the values 6 
or 7 as an upper limit for x. The importance of this 
will be seen below, since it will be shown that for n > 2 
the critical flutter Mach number decreases as  in- 
creases. 

For other values of //R and K//, this minimum may 
occur at other values of 7, but this case will not be con- 
sidered here. We shall consider the cases nm = 0 and 
n=2-—6. 

The relationship of the present work to previous 
treatments can now be observed. The major improve- 
ment over the work of Miles and Leonard and Hedge- 
peth is in the use of a shell of finite length where the 
axial modes of vibration are dictated by the end con- 
straints, and in the use of stationary rather than 
traveling waves. The end conditions actually applied 
below are those used for the beam equations, and fur- 
ther investigation is required to determine how close 
these conditions are to reality, and also to decide 
what end conditions should be applied, say, to a rocket 
structure in flight. The aerodynamic loading term is 
based on linearized supersonic theory and a static 
approach is used. The static approach is felt to be 
well justified, but while the linearized theory is cer- 
tainly valid over a wide range of Mach numbers of 
interest, if information is needed in the hypersonic 
range, the piston theory described by Ashley and Zar- 
tarian'’ should probably be used. Another develop- 
ment of the present work would be to consider per- 
turbation of a nonuniform rather than a uniform flow. 


(3) Solution of Equations 


In this section we shall formulate the complete equa- 
tions governing the flutter problem, together with the 
boundary conditions, and discuss three possible approxi- 
mate methods of attack. In Section (4) we consider 
the most fruitful of these methods in more detail and 
give an explicit solution. Since the equations for n = 0 
are basically different from those for n > 2 the two 
cases are considered separately. 

We wish to solve the equation 


: l 
L{w| = H w'(z)V, | (x —32) - Iz (3.1 
tw] J, @) (( aa) BH 


subject to end conditions corresponding to a clamped, 
simply supported tube: 


II 


w(x) = w(x) =0 at x =O) (3.2) 
P 0.4 
w(x) = w(x) =0 at x= If 
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Here L[w] is a differential operator of the fourth order 
given by Eq. (2.19) for m = 0 or Eq. (2.20) for n > 2. 
H is a constant depending on the operator chosen, 
The boundary conditions impose a condition on the 
coefficients of L, and from this condition we deter- 
mine which values of the parameters of the problem 
will lead to a complex value for the frequency. This 
will imply a disturbance growing with time—i.e., 
flutter—and the change from real to complex frequency 
will give us a definite boundary at which flutter may 
occur, 

To determine the eigenvalues of this fourth-order 
integro-differential equation in its present form would 
be very difficult. To simplify the problem the func- 
tions V,(x) are represented approximately. Now from 
Eq. (2.17) 


2m , ix 
Qi 


yo) = 
“ (a;° + n°) 


e—?t l 
(—1)” [ — 1 (3.3 
Jo p* [K,’(p)|? + [I,'(p) 2“? 


and the function occurs in the equations in the form 


d [wor (w P a) & 
dx Jo laid ‘vada BR) “~ 


where 2m is the number of zeros of K,,’(z) = O and is 
equal to 2 for nm = 2, equal to 4 for m = 3 or 4, and so 
on. If the kernel of an integral equation is in ex- 
ponential form, the equation may be reduced to a 
differential equation quite readily. Further, for n > 
0, V,(x) may be approximated quite well by the series 
of exponentials in Eq. (3.3). Therefore, if we make 
this approximation, we can reduce Eq. (3.1) to a differ- 
ential equation—but at the price of increasing the 


order by 2m. Thus, for nm = 2 we get a sixth-order 
ordinary differential equation, for m = 4 an eighth- 
order equation, and so on. 

This simplification was carried out for n = 2 and led 


to the following equation: 


d®w Qa + ae ( .) dw a a) Q2 (.) d‘*w ‘ 
dx B R/ dx° B? \R/ dx! 
L\3 dw l\4 d’w 
HIB G — PH ~ 
. (;.) dxt * | (5) dx? 


a, + a» 1\> dw , A a2 (;.) " 
G -G w= 0 (3.4) 
B (;.) dx " B YE 


The equation is solved subject to the end conditions 
(3.2) and, in addition, (since we now have a sixth-order 
equation), two conditions obtained from the integro- 
differential equation during the derivation of Eq. 
(3.4)—namely, 


d‘w dw 4 BH () dw — i 
= = at x= 
dx4 dx R}] dx? 


If a, and a are the roots of K.’(z) = 0, 


H = (16pU?/EB?)-(R/h) 
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P = amas} (1/(a:? + 4)] + [1/(a2? + 4)]} 
G = }[12/(1 — o”)j(h/R)? — 16U*(p,/E)v 


ot 
j 

We try solutions of the form w = e and this leads 
to a characteristic equation for k of the sixth degree, 


with w given by 


w C,e** + C.e** + C3e** a 
Cie** + Cse*# + Cee** 


Applying the boundary condition gives six linear equa- 
tions in the constants C;, and the condition for a non- 
trivial solution of these equations is that the sixth- 
order determinant D of the coefficients should vanish. 
This last equation, D = 0, would then give us an equa- 
tion D(/, R, h, M, v) 0 and the problem reduces to 
finding those values of /, Kk, h, and M for which the 
roots for v begin to be complex. Since the character- 
istic polynomial corresponding to Eq. (3.4) cannot 
be solved explicitly for the roots k;, the following elabo- 
rate numerical procedure would be required. 

Choose fixed values for //R, R/h, and M. Choose 
vy and find the roots k; numerically, then determine if 
D = 0. If not, vary v and repeat the procedure, until 
the two lowest roots of v are found corresponding to the 
given values of //R, R/h, and MJ. For the same values 
of / R and R/h, M is then increased and the process 
repeated until the two lowest roots for v coincide and 
subsequently become complex. This will give us the 
value of M., the critical flutter Mach number, corre- 
sponding to the given values of //R and K/h, and, of 
course, the given boundary conditions. The variation 
of v with VV is of the form shown in Fig. 7, where, in 
general, there are infinitely many values of the fre- 
quency possible for each value of \/. At some stage, 
for \f = M., two of the real values of v coincide, and 
for J > M, become complex, leading to the possibility 
of flutter. 

This method is tedious, expensive, and gets increas- 
ingly expensive as 7 increases. However, the approxi- 
mations inherent in the method can be readily assessed 
and it may be of value in testing for flutter when the 
errors in the methods to be detailed below become too 
large. An attempt was made to fit an exponential 
curve to V» in order to apply this same method for 
n 0, but a good fit was not obtainable, and the 
method of linear approximation for Vo(x) given below 
was found to be good enough for all values of the 
parameters of interest. 

If, instead of approximating V,(x) we approximate 
w(x) by means of a Galerkin or Rayleigh-Ritz procedure 

i.e., assume that w(x) is formed by superposing a few 
modes of the free vibration solution—we would still be 
faced with considerable numerical work to evaluate 
the terms on the right of Eq. (3.1). If we used the 
Galerkin approach, together with the approximation 
for V,(x) outlined above, we could proceed analytically, 
but somewhat clumsily, especially for large n. How- 
ever, we shall see in the next paragraph that another, 
linear, approximation can be made for V,,(x) that 
enables us to deal with a fourth-order differential 
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Fic 7. Values of G satisfying Eq. (4.9) for given H 


equation, and essentially to solve the problem analyti- 
cally. For this reason we shall not discuss the Galerkin 
approach further. 

Examination of Fig. 5 shows that the graphs of 
V,,(x), for moderate values of x, may be approximated 
by straight lines to quite a high degree of accuracy. 

The range of x for which the approximation is valid 
can be discussed later, but first let us consider the 
simplification it introduces into the original Eq. (3.1) 
The expression on the left is an ordinary differential 
expression of the fourth order, and is not altered by 
our approximation. For the expression on the right, 


we obtain 


d : } : l 
w'(x)V, («x — 2): dz 
dx /J0 BR 


d { w'(s , l . o { Is 
im "Cc  =_— 


where —m is the gradient of the straight line used as 
the approximating curve. This expression is 


: l 
J. BR 


l 
BR 


w'(x) — m w(x) 


where we have taken w(0)) = 0. The original integro- 
differential equation is now reduced to a fourth-order 
ordinary differential equation, and in this form can be 
solved with a minimum of numerical work. The differ- 
ential operator [L |w differs in the two cases n = 0,n > 2, 
but the essential form of the solution is the same in 


both cases. 


(4) Solution With Linear Approximation for |’, («) 


We need the following results concerning the linear 


approximation to V,(z). 


\ 


Volz) ~1 — 0.3382 OF 2 < 2.0 
Vi(z) ~ 1 — 1.62 0< 2 < 0.85 (4.1) 
Vi(z) ~ 1 — 2.52 0<2< 0.65 


Since each function V,(z) appears in the equation in 
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the form V,[(x — z)-(//BR)], with x — z ranging 
from 0 to 1, the approximations will be valid provided 
that //BR is less than 2.0, 0.85, or 0.65 for n = 0, 4, 
or 6, respectively. We shall use the relation 

V,(z) = 1 — mz (4.2) 


and substitute for m at a later stage. 

For n > 2 we assume a solution for w in the form 
w = w(x)-cos n6-e'”. Substituting in Eq. (2.20) and 
using Eqs. (2.18) and (4.2) we obtain 


d*w a -__ alt Iv Oe . 
det + (;) d?(n? — 1)?n4w — n‘ (5) E w= 


(5) R? pU* {dw l 
nN : . q ‘ WwW 
R/ Eh Bl \dx BRS 


This can be rearranged to give 
(d4w/dx*) — H(dw/dx) — Gw = 0 (4.3) 


where H = (n‘pU?/EB)-(R/h)-(l/R)* 


and 


— «py? — dn? — 1)? — 


mpU? (*) 
EB? h 


We shall use the end conditions for a rod clamped at 
x = 0 and simply supported at x = 1 


then w= dw dx ins 0 at x=0 | (4.4) 
w= d*w/dx? = 0 at x= 1 f 


We wish to determine when v becomes complex for 
varying values of R/h, //R, and M (or B). This is 
equivalent to finding when G becomes complex and we 
regard G as an eigenvalue. At the end of the solution 
we can determine v from the resulting value of CG. 

We assume a solution of the form e“” in Eq. (4.3) and 
obtain the characteristic equation 


kt— Hk-G=0 (4.5) 


This, in general, has four distinct roots k; (¢ = 1, 2, 3, 
4) so that the general solution may be written in the 


form 
w(x) = Cal + Col + Col + Ce” (4.6) 


If the values of G and H are such that some of the roots 
coincide, we need only vary one of the parameters 
R/h, l/R, B, or v slightly, to arrive again at the general 
case of four distinct roots. We shall therefore assume, 
without loss of generality, that the solution may be 
written in the form of Eq. (4.6), with all the values 
of k distinct. 

We now impose the boundary conditions (4.4) and 
derive four linear equations for C;, the coefficients of 
which are functions of G and H. The condition for a 
nontrivial solution is that the determinant formed by 
these coefficients should vanish, i.e., D = 0, and this 
equation then gives the relation between the param- 


eters of the shell and incident field of flow to deter- 
mine under what conditions the values of v become 
complex and flutter develops. 

To solve the characteristic equation (4.5) explicitly 
for k; in terms of H and G is of course possible, but 
inconvenient, and another method is used here, follow- 
ing Hedgepeth and Stepanov. In this method, instead 
of solving for all k; in terms of H and G, we essentially 
solve for 7 and three of the roots in terms of G and the 
fourth root. From the symmetric properties of the 
roots of Eq. (4.5) it is clear that these may be written 


in the form 


be otf h=- a-#\l as 
ky = -—a+t+ ly ky = —a — iyf a 


The sum of these roots is zero, as required by the form 
of Eq. (4.5). 


In addition we have the following results: 


7 — & — 2a* = 0 
2a(y? + 6B?) = H 
(a? — B?)(a? + y’?) = -—G 


Solving for 8, y, and G in terms of a and H we find 


| (= ‘) 
g2 = es gail 
4a 
; H a 
7 = ( + a’) (4.8) 


— 4a‘ 


16a? 


When we form the determinant corresponding to the 
boundary conditions (4.4) we obtain 


l l | ] 
re 2 & & law 
é e 


e* e 
ke" ko2e*? ks2e"* hye" | 
Substituting for the k; from Eqs. (4.7) we get 


2a8y sinh 2a + B(8? — a?) sin y cosh B — 
v(y? + a’) cosy sinhB = 0 (4.9) 


Due to the nature of the equations we are unable to 
obtain an explicit algebraic equation for G (or v) and 
it is necessary to solve this transcendental equation 
numerically. The geometry of the cylinder is kept 
fixed and H/ and G are used as variables rather than / 
and v. We therefore fix H, and vary a until Eq. 
(4.9) is satisfied—this gives a value of G corresponding 
to a particular value of 7. This process is repeated 
until the two smallest (positive) roots of G are obtained. 
We next choose a larger value of HT] and repeat the 
process until the roots of G coincide at HW = Heit. 
From this value we then determine /,,;, for the given 
geometry of the cylinder. 

In practice we do not actually evaluate the roots of 
G at each stage, since we are only interested in the 
value of H giving coincidence. Hence for each value 
of 77, sufficient values are taken to ensure that two 
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distinct roots still exist. As the value H = H.,j is 
approached, we take smaller steps for H and a, and 
seek a more exact result. The results are given in Fig. 
7 and will be discussed in the next section. 

For n = O we have a similar equation to solve. 
From Eq. (2.19), using Eqs. (2.18) and (4.2) we get 


dw (;) 1 l4p,hU? oe 
ix! R) a — 
lt pU? {dw l 
. —m We 
D Bi \dx BR S$ 
or (d‘w/dx*) — H'(dw/dx) — G'w = 0 (4.10) 
where 


H’ = 


12(1 — 6*)pU? a] 
EB h 


: (:) [eee os ee) 
G = ~— —m . 
h D ad? B*D 
9 


The solution is exactly the same as in the case n > 2, 
but the coefficients H’ and G’ have different values, 
and this will lead to different critical Mach numbers 
for nm = 0, as might be expected. 

The effect of axial and circumferential tensions can 


readily be obtained for the case m = 0. When N, 
and N, are not zero, Eq. (2.19) is replaced by 
o'w 4U*p.h Ow Ph Ow 


ox! RD a  D "dx 


(;.) l I 
-_—WY = Z 
R d* D 


and after some reduction, this gives, in place of Eq. 
(4.10) 
d*w 


dx! 


F = (I?h/D)-N, 


where 
This equation can be handled in much the same way 
as that applied to Eq. (4.3) above. Again, we may 
assume roots of the form 


ky» = a’ ~ B’ ks: = —a’ — iy’ 


but now, instead of Eq. (4.8), we have 


: Hy’ ; F 
ss 6S ~— gg’? + 
ta 2 


z H’ | F 
lied eae 


| 
G’' = 

16a’? 
Since the boundary conditions are the same as before, 
the determinant D is unchanged, and the final equa- 
tion is still Eq. (4.9), with primed quantities intro- 
duced. As F increases, /7’.,;, increases, and if F 
becomes negative so that the cylinder is under axial 
compression, then //’,,;, decreases. 


Pi 


\ 
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(5) Results and Discussion 


To provide a basis for comparison, the results are 
first drawn in Fig. 3 using the same values of param- 
eters as those quoted by Stepanov. Later, in Fig. 8, 
the results are shown for a range of parameters that 
may be more realistic, in terms of the equations used 
here. Fig. 3 shows the curve of critical Mach number 
against thickness ratio for a steel cylinder at sea level, 
with the end x = 0 clamped, and the other end simply 
For comparison with Stepanov’s results 
we have taken giving four circumferential 
waves, and the length ratio //R 10 and 12. 
readily be seen from the diagram, the thickness ratio 
required to prevent flutter at a given speed is somewhat 
smaller for the results based on our new aerodynamic 
term, and this result holds for all comparisons between 
our results and those given by Stepanov. In other 
words, the more realistic aerodynamic terms have ex- 
tended the flutter boundary. as mentioned 
earlier, it is difficult to know which value of is really 
significant, and the result shown in Fig. 3 should not 
be interpreted as giving absolute figures to prevent 
flutter at all modes, but only to be used if the circum- 
= 4. 


supported. 
n= 4, 
As can 


However, 


ferential modes can be so constrained that 


It is interesting to compare these results with those 
obtained by Leonard and Hedgepeth for the case of an 
infinite cylinder. The Leonard and Hedgepeth curve 
for the case n = 4 (or 5) is shown in Fig. 2, and the 
significant effect of finite length is at once evident, 
since, for the same speed, the thickness ratio is now 
only about one tenth of the value given by Leonard 
and Hedgepeth to ensure stability. 


One final comment is relevant before proceeding to 
a more detailed analysis of the results obtained in this 
paper. For the infinite cylinder there is no significant 
difference in the answer given in Fig. 2 from that of 
Miles, with n = 0, and Leonard and Hedgepeth, with 
n = 4or 5. For the finite cylinder, as we shall see, 
this agreement no longer exists and there is a large 
divergence between the curves for » = () and those 
for larger values of n. This again shows the difficulty 
in extrapolating results obtained for vibrations of in- 
finite cylinders to the case of finite cylinders, with end 
conditions playing a vital role, both in the formula- 
tion of the problem and in the results obtained. 

In view of the remarks made earlier concerning a 
possible upper limit for ”, we have plotted the results 
for nm = 6 in Fig. 8, again for a steel cylinder at sea 
level. with //RK = 4or6 

The equation used to determine the relationship 


between J. and h/R is 


n'pU? (*\(2) 5 
EB \h/\R/ 


and for the end conditions clamped and simply sup- 
Using the equation 


H = (5.1) 


ported, Fig. 8 gives [7/.,;, = 480. 
of state for air we can put this in the form 
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Fic. 8. Variation of M, with thickness ratio for a steel cylin- 
der at sea level, m = 7. End conditions: clamped—simply 


supported. 
niyp ( | \G) (;.) 
B . = 480 
E ¥ B/\h R 


where p = air pressure in the free stream. Separating 
the velocity terms on the right we obtain 


(2 = 5) = 480 < (=) () (5.2) 
BJ) niyp \R l ies 


In numerical work we have taken 

| = 2 X 10! dynes/cm.’ 

OY = 1.4 

Psea-tervet = 1.014 XK 10% dynes ‘cm.’ 
We shall analyze each of the terms in Eq. (5.2) sepa- 
rately, to determine their effect on the flutter boundary. 

The speed or Mach number only occurs in the form 
(B + 1/B) where B = (M2? — 1)'”*. The term (B + 
1/B) has a minimum value equal to 2 for B = 1 (i.e., 
uM = V2), and this implies that we are unable to in- 
vestigate the region | < M < V2 with the present 
theory. If any such attempt is made, it gives the 
anomalous result that, for a certain value of thickness 
ratio, the panel would be unstable between J/ = | 
and M = M, < V2, then stable until the speed 
reached the second, real, critical Mach number > 
V2. This restriction is essentially forced on us by 
using linearized supersonic theory, which o/ course 
breaks down in the lower supersonic region. In our 
case this region extends to M = V2. Itis perhaps of 
interest to note that Leonard and Hedgepeth found 
that the use of an infinitesimal amount of structural 
damping for the infinite cylinder apparently removed 
this restriction from their results, and enabled them 
to proceed to M = 1, even though they used linearized 
supersonic theory for the aerodynamic loading. The 
difficulty in the lower supersonic region is probably a 
consequence of using the static approximation and 
could be avoided if time-dependent terms were taken 
fully into account. 
Although no qualitative results were developed in 

the case of structural damping, an analysis of the 


equation showed that there would be no discontinuous 
change in the results if damping had been introduced. 

Further, unless we are particularly concerned with a 
vehicle traveling for extended periods in the range 
L< We< V2, there seems little merit in developing a 
theory to investigate this range, since in practice most 
rocket shells will pass through this range in such a 
short tine that panel flutter is unlikely to have any 
significant effect. 

The significance of the particular circumferential 
mode of vibration is seen to be very great, since it 
enters as a term ~*. We have already discussed 
the range of possible values of m, and we may conclude 
that this is a sphere in which much more work, both 
experimental and theoretical, needs to be done. 

The free-stream pressure appears in the denominator, 
so that a decrease in pressure will lead to an increase 
in \/, for the same thickness ratio. The free-stream 
pressure will depend on the total configuration of the 
body and hence it is easier to make comparative rather 
than absolute statements. Since 


Psea-tevet ae 10 P50; 000 ft = 100 P100,000 ft 


it can be seen that height also has a very important 
effect on the flutter boundary. (See Figs. 9 and 10.) 
Again, pressure plays a very much greater role in the 
equations for a finite cylinder (both in Stepanov’s 
results and the present ones) than for an infinite cylin- 
der, where Leonard and Hedgepeth and Miles find 
that the speed of sound, rather than air pressure, is the 
decisive parameter representing the effect of height. 

As might be expected, the length ratio also plays a 
dominant role in the equation, for \/, is approximately 
proportional to (X//)*. The thickness ratio only ap- 
pears to the first power, but of course, in design and 
fabrication of these shells, this ratio is probably the 
most important. 

The material of the shell again plays a more impor- 
tant role than in the infinite cylinder case, where its 
only effect is in the term (1 — o*), where o is Poisson's 
ratio, and does not vary a great deal in diferent mate- 
rials. In our case, the properties of the material enter 
the flutter equation as X, Young's modulus, which can 
change appreciably from one material to another 
e.g., Eai:o (1/3) E steer. 

For the other case considered, n = 0), the equation 
obtained to determine the flutter boundary, analogous 


to Eq. (5.1) is 


fy = me #3 U? (1) = 480 (5.3) 
~  -«—EB wn atin | 


After some reduction, this gives 


(2 é l ) 480E (| ) (7) - 
= . . (5.4) 
B 12.1 — o*)yp \R l 


This equation is similar to the equation for n 2 = 
[Eq. (5.2)] in the dependence of 1/, on Young’s modu- 
lus, the air pressure, and the inverse cube of the length 


ratio. Now, however, the dependence is on the cube 
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of the thickness ratio, rather than merely the first 
power. 

The approximation chosen for V, earlier [Eq. (4.1) ] 
does not appear to have a critical influence on the flutter 
boundary. In fact, the slope of the straight line ap- 
proximating V’, does not enter the expression for the 
critical Mach number at all, and the only way I’, in- 
fluences this equation is by its initial value—which 
in each case is the same—unity. The term involving 
the slope (—m) does affect the flutter frequency in the 


term 


.. . FL] Ops’? oar . mpU?/R 
G=n*ti— — d*(n* — 1)? — 
R E EB? \h 


It will be recalled that the two smallest roots of G 
(or v) coincide at 7 = H,,;., and by equating G to the 
numerical value of this coincident root, we can get the 
frequency of the cylinder at the onset of flutter. Thus, 
if we let this value of G be G,,;;, we have 


E J 1 /R 


4 
oe (7) Gorie + d2(n? — 1)? + 


mpU* (7) a 
- rp (5.5) 
EB? \h/S 


for the case n > 2, and 


D iG) @ 4 l 4 mpU*R*) 5.6) 

y? = rae - (5.6 
psUhR? \\ 1 '" @ BD J 

for the case nm = 0. Hence the flutter frequency may 


be obtained, if desired, in any particular case. 

At this stage we can investigate the validity of the 
range of approximation for V(x) given by Eqs. (4.1). 
For n = 0, //BR must be less than 2.0. If //R = 4, 
this implies that B > 2 or M > 2.2, and as (//R) in- 
creases our results strictly apply only to a region be- 
ginning with an increasingly high Mach number. For 
n = 6, the approximation is good for M > Gif (//R) = 
4. However, only the initial value of V(x) influences 
the flutter speed, since on physical grounds the func- 
tion V,(x) acts to damp out the effect of the aero- 
dynamic load [see Fig. 5 and Eq. (3.1)|, we feel that 


Mc 


Variation of J/,. with thickness ratio based on Eq. 
(4.10), for a steel cylinder at sea level, n = 0. 


Fic. 9. 


ANEL FLUTTER 207 


& 
o2 
e 
78 
—_— ; 
£ 
© = a 
R _ 
F 2 3 « 5 Mc 
Fic. 10. Variation of M/, with thickness ratio based on Eq. (4.10) 


for a steel cylinder at 50,000 ft., » 0 


the results will not change discontinuously outside 
the range given above. Further, for all reasonable 
values of M/ (say, <30), it is seen by comparing Figs. 
S and 9 that ” = 
thickness ratio required to prevent flutter at a given 
Mach number must be considerably greater for 7 0 
than for the case » = 6. We are therefore more inter- 
(), and it is here that the range of 
It is of interest to note that the 


() is the more critical case—i.e., the 


ested in the case 1 
validity is greatest. 
wavelength in this critical case is slightly smaller than 
the cylinder length. 

All the previous discussion refers to the particular 
end conditions used in our example, but it is important 
to note that it applies equally well to any set of end 
conditions. The effect of varying the end conditions 
is to change the numbers in Fig. 8, but not the form of 
the curves. Instead of H,,; 180 obtained in the 
above case, the case of a cylinder clamped at both 


ends lead approximately to a value /7/,,;; (3/2)-480, 
and for the case of a cylinder simply supported at both 
ends Hw. = (3/4)-4S0. Results for these cases 


could then be obtained easily. 
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Nonequilibrium Flows (Continued from page 196) 





This next approximation has not yet been worked 
out. The following is an outline of how this might 
be done. By developing gradient functions for equilib- 
rium flow, the curvature of the equilibrium shock wave 

Then let the actual shock have this 
For this curved shock, a_ shock- 


could be obtained. 
same curvature. 
oriented coordinate system could be used to attempt 
to get the next approximation to the transition region. 
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Dynamic Airloads and Aeroelastic 
Problems at Entry Mach Numbers’ 


GARABED ZARTARIAN,* PAO TAN HSU,* ann HOLT ASHLEY* 
Massachusetts Institute of Technology 


Summary 


During the process of entering planetary atmospheres, space 
vehicles may encounter their largest flight dynamic pressures at 
very high Mach numbers. The most severe aeroelastic phe- 
nomena may therefore occur at hypersonic speeds, and a need 
exists for reliable methods of dynamic airload prediction in this 
range. Two approximate techniques are described, which hold 
promise of supplying such information for two- and_three- 
dimensional pointed bodies of rather general shape. 

For low to moderate values of the parameter /,,7, a vari- 
ational principle is used in conjunction with hypersonic small- 
disturbance theory (thickness ratio r<_ 1) to find the pressure 
distribution due to a steady or arbitrarily time-dependent motion. 
The key step is a reduction to a series of equivalent two-dimen- 
sional piston problems, which are solved by the Ritz method. 
This procedure has been successfully applied to wings and sur- 
faces of revolution, and the generalization to more complicated 
shapes is discussed. 

At higher values of J/,,7, a simple extension of shock-expansion 
theory, also relying on the piston analogy, appears to furnish a 
rational approximation for many configurations. It consists of 
determining flow parameters just aft of the nose shocks, then 
carrying out a Prandtl-Meyer expansion along surface stream- 
lines in planes normal to the body. For thickness ratios typical 
of high-L/D entry vehicles, dynamic loading can be calculated 
up to any Mach number of aeroelastic importance provided the 
leading edges are not too blunt. 

Several illustrations of both theories are given, including com- 
parisons with more exact aerodynamic results and application 


toa case of chordwise divergence. 


Principal Symbols 


speed of sound 


a 

b = semichord of airfoil 

i. 1,% unit vectors in the x-, y-, 3-directions, 
respectively 

k = reduced frequency, k = whb/l 

l = body length or characteristic length 

m = rate of change of Mach number with re- 
spect to body inclination, (just outside 
the vortical layer for three-dimensional 
bodies ) 

m = rate of change of Mach number with 
respect to body inclination under the 
vortical layer for three-dimensional - 
bodies 

p = pressure 
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q = dynamic pressure 

q = velocity vector 

,, H, 2 = cylindrical polar coordinates (or r, y, ¢ 

t = time 

to = time when a given fluid slab enters the 
apex of the configuration. Identifies 
fluid slabs 

hi, te = starting and terminal times for a given 
fluid slab as it travels over the body 

u,v, W = components of fluid velocity in x-, y-, 
directions, respectively 

Up piston velocity 

Ur = radial velocity 

X,Y, 2 = Cartesian coordirates 

Sis = piston location (one-dimensional cases) 

Zz = unsteady motion in the z-direction 

21 = amplitude of unsteady motion z = 2,f(t) 

a. on = unknown coefficients in assumed series for 
the velocity potential 

A.C = aerodynamic center for airfoils 

B(x, y, z, t) = 0 = equation of body surface 

cs = pressure coefficient: 


Cy = (~ — fe)/(1/2pq,U*) 


I = variational integral 


L = lift (positive up) 

M = Mach number 

R or R, = radius of body or piston 

S(x, y, 3, 1) = 0 = equation of shock surface 

U. = free-stream velocity 

V = disturbed volume of fluid 

a = angle of attack; also pitching motion with 
amplitude @ (positive nose up) 

ao = steady mean angle of attack 

B = shock-wave angle 

7 = adiabatic gas constant = C,/C,; = 1.4 un 
less specified otherwise 

p = density 

T = thickness ratio 

9 = velocity potential 

w = circular frequency in radians per second 

9 = bar indicates ‘‘unyawed”’ quantities 

( dn = subscript N indicates quantity at nose 

( Yor Ja = subscripts U, L indicate upper and lower 
surfaces, respectively 

_ = subscript © indicates quantity in the 


freestream 


Introduction 


. IS WELL-ESTABLISHED that nonlinear aerodynamic 
effects, such as those of wing thickness and body 
cross-sectional shape, can significantly alter aeroelastic 
stability even in the lower supersonic range. Pro- 
gressive increases in the hypersonic parameter M.7r 
will intensify these nonlinearities. Experience also 
suggests that the most severe aeroelastic problems occur 
within the region of highest dynamic pressure g., along 
the flight envelope of a hypervelocity vehicle. (To be 
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more precise, for an envelope or flight path which con- 
tains a constant-q., segment well above transonic speeds, 
the trouble spot usually lies at the low-Mach-number 
end of this segment.) As can be deduced from many 
examples in Chapman’s paper,! typical entry tra- 
jectories display the high-g., point around J/., = 10. 
Similar peaks in the hypersonic range can be anticipated 
for rapidly-accelerated vehicles on the way out of the 
atmosphere, for instance, the antimissile missile. 

The need cannot be questioned, therefore, for a 
manageable hypersonic airload theory comprehending 
rather general shapes of wings and bodies, which deform 
in an arbitrary steady or time-dependent fashion in 
response to the loads. Concerning the case of unswept 
and moderately swept airfoils, considerable effort has 
been expended to estimate the influence of profile shape 
in two-dimensional unsteady flow (references 2 through 
6, among many others). For slender bodies, one ex- 
pects linearized methods to remain valid up to higher 
values of \/..7 than for wings. Moreover, even simple 
Munk-Jones theory accounts for body thickness 
through the dependence of crossflow virtual mass on the 
geometry of sections normal to the flight direction. 
These observations are indeed fortunate, because the 
few solutions of the three-dimensional nonlinear field 
equations (cf. references 7-12, 15) are generally limited 
to steady flow over specific configurations, and the in- 
fluence of a blunt nose on the pressure distribution aft 
is only poorly understood. 

Two exceptions to the steady-flow limitation occur 
in the works of Eggers and Savin,!* who indicate the 
possibility of adapting their shock-expansion technique 
to oscillatory motion by following particle trajectories, 
and of Holt.'* Some discussions of unsteady hyper- 
sonic flow are also given in the book by Hayes and 
Probstein,’ in particular, a thorough examination of 
the shock layer on a blunt surface. The present paper 
emphasizes configurations whose surfaces have small 
slopes almost everywhere. It is, in a sense, comple- 
mentary to treatments of the shock layer, and there 
exists an important, unfulfilled need for means of 
fitting together the two types of approximate solution. 

In what follows, there are described two techniques 
that hold promise of supplying the required hyper- 
sonic airload information for a variety of pointed or 
slightly-blunted shapes. An inviscid fluid is assumed, 
which will cause appreciable errors at high J/..’s and 
low densities but is subject to adjustment by adding a 
boundary-layer displacement thickness to the body 
surface. In a similar way one can conceive of a thin 
vortical layer of high-entropy fluid, originating from 
the detached bow shock and spreading out over the 
surface. This mechanism offers a hope of adapting 
slender-body theory to account for slight nose and 
leading-edge bluntness. It should be more successful 
in the case of three-dimensional vehicles, but much re- 
mains to be investigated before practical applications 
will be possible in unsteady flow. 

The starting point for both of the present techniques 
is Van Dyke’s extension'® to pointed bodies of Hayes’ 
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hypersonic similitude.” The fluid motion is shown to 
occur principally in planes normal to the flight direc- 
tion, so that wing and body problems are equivalent to 
one- and two-dimensional pistons pushing into gas 
which was initially at rest. The easier method to apply 
is developed first. It is a simple variant on well-known 
shock-expansion theory, whereby flow quantities are 
determined just aft of the shock attached to the pointed 
nose or leading edge and Prandtl-Meyer expansions 
then carried out along predetermined lines at the body 
surface. For two-dimensional shapes, the entire range 
from moderate to very high /., can be covered. For 
figures of revolution and the like, however, the pro- 
cedure proves to be justifiable only when M/., > | 
and the product J/..7 exceeds a minimum value near 
unity (or higher, depending on body shape). Thus, a 
gap remains, in the case of three-dimensional slender 
configurations, at the lower end of the W/.r-scale. 
This is filled by the second technique, which consists of 
solving, by the Ritz method, a variational statement 
of the analogous two-dimensional piston problem. 
Isentropic flow is assumed, thus putting firm upper 
limits on the permissible J/..7, but it is believed that 
this restriction can be removed by recourse to more 
general variational principles (e.g., Skobelkin'’’ and 
Fto**).. 

A word of justification is in order for going to the 
variational approach. It has, of course, been used 
with great success in the field of elasticity, but has 
seldom found favor in fluid dynamics, where compli- 
cated outer boundary conditions and tedious computa- 
tions are encountered. The works of Wang, ef al.'9~*° 
for subsonic steady flow, Lieber and Wan,*** Fyfe 
and Klotter,*6 constitute most of the published applica- 
tions. Among these, only one* deals with unsteady 
flow, specifically with oscillations in a one-dimensional 
gas-filled tube. Their paper recognizes the advantages 
of the method that recommended it for the present 
analysis: the ease with which nonlinear phenomena are 
described, the particular suitability for fields which are 
finite in extent and free from singularities. No mathe- 
matical proofs of convergence or uniqueness have been 
found. Therefore, a number of examples are given to 
establish the validity in practice. These also serve to 
clarify the main features of the rather complex pro- 
cedure, and to guide the reader toward the ultimate 
goal of treating quite arbitrary unsteady boundary 
conditions and sectional shapes 


The Piston Analogy and Equivalent Body Shapes 
for Unsteady Flow 


The hypersonic small disturbance theory of Van 
Dyke" is based upon an examination of orders of mag- 
nitude of the terms in the field equations, surface 
boundary conditions, and shock-jump conditions de- 
scribing flow past a slender, pointed configuration in a 
stream of perfect gas at M.. > 1. Without going into 
the full mathematical details, consider small, time- 
dependent lateral motions of the system, shown in 
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Fig. 1, which is described in inertial coordinates whose 
x-axis is aligned with the free stream and (as close as 
possible) with the body axis. A boundary condition is 
applied on the normal component of fluid velocity at 


the surface 
B(x, y, 2, t) = 0 (1) 


The familiar oblique shock equations govern the dis- 
continuities in velocity q and properties p, p, s across 
one or more shock surfaces 


S(x, y, 3, t) = 0 (2) 


The field differential equations require conservation of 
mass, momentum and particle entropy between these 
bounding surfaces. 

Reference 16 introduces dimensionless independent 
and dependent variables, which are of order unity 
throughout the disturbed flow, by means of the trans- 


formations 


Z=x/l, J = (1/r)(Cv/), 
= (i/ris/), t= US 


(3a-d) 


and 


= 7 TD, 


[1+ 774], v= U.77, w= 


U 
p=p [yM 277] D, p= p.p (4a-e) 


Here / is a streamwise characteristic length, and 7 is the 
thickness ratio, assumed to be of the same order as the 
maximum body slope. The requirement that p = 0 (1) 
sets a lower limit on the permissible size of M.r 
for three-dimensional shapes like the one in Fig. 1. 
Van Dyke’s development can also be adapted" to 
wing-like configurations by dropping dependence on the 
spanwise variable y, but, in such cases, it can be shown 
that \/..r may become arbitrarily small, provided \/ 
is large enough. 

When the transformations Eqs. (3)—(4) are employed 
and terms 0(r*) are dropped from the hypersonic-flow 
equations, an important simplification occurs: the 
streamwise velocity @ is decoupled from all the other 
unknowns, and the remaining equations contain the 
quantities @ and 7 only as the operator [(0/0/) + 
(0/0%)|. Thus the number of independent variables 
is reducible by one through the transformation to still- 


air coordinates 


J =7-i, tP =i V=V7, 2 =2 (5a-d) 


A new, equivalent unsteady flow problem emerges 
which consists of two-dimensional motion in the 7’ — 2’- 
plane produced by an expanding and contracting im- 
permeable boundary. Its equation, obtained by sub- 


stituting Eqs. (3) and (5) into Eq. (1), may be written 
B'(y’, 2’, ’) = 0 (6) 


Physically, the problem is that of the transient dis- 
turbance created by a two-dimensional piston, whose 
motion commences at some initial instant 7’ = /)’ when 
the particular fluid slab strikes the body’s nose. If 
the original motion is itself unsteady, so that B de- 
pends explicitly on /, the piston behavior is different 
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Fic. 1. A slender configuration in hypersonic flow 


for each slab, and B’ in Eq. (6) depends parametrically 
on fo’. Hence, a number of piston problems has to be 
solved in order to compute a single case of transient or 
oscillatory loading on a hypersonic vehicle, whereas 
only one is needed for a steady flow. 

In the shock-expansion theory which follows, it is 
convenient to employ the piston analogy in reverse. 
That is, the expanding-contracting piston in any one 
slab has as its counterpart the steady motion of a dis- 
torted body 


B’( 5’, 3’, '; to’) > Bz, 5, 3, lo + FB) 


x (7) 
B (« y, 8, to + U ) 


For example, consider a body of revolution with radius 
R(x), whose axis performs small lateral excursions or 
bending deformations described by 


3 = 2(X, bt) S) 
Eq. (1) for this case reads 
B(x, y, 2, t) = R*(x) — iy? + [s — a(x, t)]?} = 0 (9) 


O att = fo, the piston 
equation (in dimensional, still-air coordinates) is 


For the slab striking the nose x = 


B'(y’, 2’, t’') = R*[U.(t' — te)| - 
iy’? + [s’ — 2,(U.(t’ — to), t’)|?} = 0 (10) 


This slab ‘‘sees’’ a body in steady flow whose equation is 


B(x, y, 2; to) = 


R*(x) — \" +ist= a(x, to + + yy QO (11) 


Eqs. (11) and (10) are used, respectively, in connec- 
tion with the shock-expansion method of the next sec- 
tion and the variational solution of the piston problem 
itself. 


An Extension of the Shock-Expansion Method 


<xtensive and successful applications have been made 
of the shock-expansion method to pointed bodies in 
steady-flow.'*: *~*° For the calculations of surface 
pressures, it consists simply of (a) determining the flow 
conditions at the nose just aft of the shock from exist- 
ing solutions for the nose tangent-cone [or the leading 
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Fic. 2. Geometry of a two-dimensional airfoil in unsteady 
motion. 


edge tangent-wedge for airfoils], and (b) using Prandtl- 
Meyer expansion from the nose along certain paths on 
the surface. Step (b) is straight forward if no signifi- 
cant compression turning is called for along the body. 
Should there exist appreciable compressive slope dis- 
continuities, the resulting shocks must be accounted for, 
and the expansion carried out stepwise. When com- 
bined with the equivalent steady-state-body concept 
of the hypersonic small-disturbance theory, it can be 
extended to cover unsteady motions at high Mach 
numbers. Being an approximate procedure, it has 
certain limitations which will be pointed out in later 
discussions. 

As an introduction to the more complex situation 
for bodies, the two-dimensional oscillating airfoil will 
be analyzed. For low values of the similarity param- 
eter \/..7, say up to 0.3, the second-order piston theory 
of Lighthill® has been employed for many aeroelastic 
applications.*': ** His third-order solution, although 
derived on the assumption of isentropy, is expected to 
hold for somewhat higher values of J/..7 despite the 
fact that appreciable entropy changes occur at the lead- 
ing-edge shock.* In what follows, the leading-edge 
shock effects are included in an approximate manner. 

Consider the thin airfoil depicted in Fig. 2, which 
performs small lateral oscillations about some mean 
angle of attack a». Let zy(x) describe the ordinate of 
the upper airfoil surface relative to the chord line and 
z,(x) that of the lower surface. Since the top and 
bottom surfaces do not interact, it is permissible to 
treat each side independently. If the motion of the 


chord line is given by 

zi(x, t) = 21(x)f(t) (12) 
the equation of the upper surface becomes (for ay < 1) 
B(x, 2, t) = 2 — [gu(x) — av + a(x, t)] = 0 (13) 


Following the discussion of the previous section, a fluid 
slab located at time ¢ = ¢) at the nose is located at time 
tat x = U.(t — to). Furthermore, this slab experi- 
ences the same disturbances as if it had passed in steady 
flow over a warped upper surface 


My a (7 + 1 
M, 


a eee 


* For the lower surface this angle is {(dzz /dx) my + ay}. 


2 {(M.8)? — [(y — 1)/ 


B(x, 2) = 
z— {ay(x) — aov + 21(x)f[to + (x/U.)]} = 0 (14) 


The boundary condition is 


q:-VB = 0=2 (U.i + wk)-VB, on B=0 (15) 


since wu = U.,[1 + 0(7r?)| & U.,; the surface flow inclina- 
tion 6(x), relative to the x-axis is (to the order of accu- 


racy of this analysis) 


ae ] a sdzy 
6(x) = U w= . a ed —_— 
fda x -.. ltr + (x/U-))\ 
ie I (« + U )+ 21(X) we ( (16) 


= {On(x)} + te(x, to)} 


which takes the value at the nose (x = 0) 


1zy 
by = 6(0) = (32) — ay : os 
aX J r=0 


iz). Of lto + (x/U.) 
J f(t.) + 2:(0) & | | 
Ox 


— 
ldx r=0 r 5 oF 


= {Ay} + {ex(to) 
Oy represents the mean inclination of the upper surface 
(at the nose) for all fluid slabs, and €y(fo) the perturba- 
tion angle which is dependent on the particular fluid 
slab. Oblique shock formulas for large Mach numbers 
and small turning angles give, for the pressure at the 


nose 
Pn/po = [27/(y + 1) [(6M.)? — 1] +1. (18) 


where @ is the shock inclination corresponding to the 
turning angle é6y. For small time-dependent parts of 


the motion ey, 8 may be expanded about its value 6 due 
to On 


8B = Bll + (1/)B(0B/08) ey] (19) 


In a similar way, one can write the pressure py and 
Mach number J/,y behind the shock in terms of the 


values Py and J/y, they would have if, ey = 0 
bu = pull + nen] (20) 
My = My[1l + mey] (21) 


Here 7 and m are the dimensionless rates of change of 
pressure and Mach number with turning angle 6. Ex- 
pressed in terms of the known angle 6y = (dzy/dx);~» — 
ao, * these various quantities can be shown to be 


9 


: r 2 
M.B = _ (M..6x) + 4 + M os | 44 
(22) 


by/b. = [2¥/(v + 1)](M.B)? -— (vy — D/(v + DI 


(M..B) , 
4 o iz ‘ 2 pee She : V1/2 (24) 
2)}"/? {[(v — 1)/2] (M8)? + 1}? 
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Pe = (7 + 1) X 
(M.B)* : 
{y(M..B)? — [(y — 1)/2]}[(7.B)? + 1] 7 
“ la , (v¥+)) 
Mo 062M 2 
M.B) 
(26) 


{ly — 1)/2] (7.8)? + 14 [(07.8)? + 1 


Once the nose conditions are determined, the pres- 
sure at any point on the surface described by Eq. (14) 


p(x) 7 Pw Pn P(x) 


PN p Py Py 
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to 
w 


can be computed by the well-known Prandtl-Meyer 


expansion formula: 
p(x) ec: Aocall . Pow 


All these expressions satisfy the hypersonic similarity 
rule. Eq. (27) is applicable only at high Mach num- 


(27) 


bers and low values of [6(*) — 6y| and is consistent 
with the approximations involved in Eqs. (22)—(26). 
Combining Eqs. (16), (17), (20), and (21) with Eq. 


(27), one obtains 


) (28) 

A, ae 27 l 

Be % t {1 + nex} }! + [(y — 1)/2]My[1 + mey| | (ac) — ox) («x to) — v(t) | 
2 

Eq. (28) may be expanded to yield, [with A(x) = 6(x) — @y and O(x, to) = €(x, to) — en(to)] 

p(x) Ps) (? — *) pa ern - men A(x) + O(x, to) } = 

= 41+ MyA(x) ¢ sl + + yMy = > + O(e?) (29) 
p (5 { 2 . f { = iis {1 + [(y — 1)/2] WyA(x)}1 ——_ | 


For small lateral motion, the terms 0(¢°) may be neglected. 


Finally, using subscripts LU’ and L to identify quantities 


on the independent upper and lower surfaces of the airfoil, one can derive the load distribution 


Pt — Pu yM 


- , (Cp, pit Cr, 





9. 


Pru J y¥-— 1 - - 7 
¥ 41+ ( 2 ) Avedetx)| 


+ [(y — 1)/2]MwrA,(x) 5 


= muyenAy(x) + O(x, to) 
1+ 7 + yMnv — 
IUEN 7 NL iy we l(y— 1) 2) M wu Ay (x) } 


Day —1\ _ si it ss myenAz(x) + O(x, to) a 
“= (‘ ait +- € “ ) MwA, (x) ;! — new — YM yx iy a — (30) 
4 - I 


For a symmetrical airfoil at zero initial angle of attack Myr = Myv, Bt = Bu, ete. Then a number of cancellations 
yields the simpler result 
Pr, — Po Px 771 » alias allie j y¥-1. a ) 
me i 1+ : M vA(x) nex 1+ : MyA(x) | + yMy|[mA(x)ex + O(x, to)] 
p b Z 2 
(31) 


In order to obtain the complete loading associated with the oscillatory motion, the process just described must be 


repeated for all fluid slabs throughout one period of oscillation. 
at any time that all points x on the surface of the actual airfoil. 
Since Eq. (31) is in closed form—.e., fo is left arbitrary, fo can be replaced 


by a series of slabs with fp = ¢ — (x/U.). 
by ¢ — (x/U.). 


In particular, it is necessary to find the pressure 
This involves the pressures experienced at time /, 


With this substitution, Eq. (31) expresses the desired pressure distribution at any time ¢. (It 


should be emphasized that this simplification is possible because of the introduction of the approximate expansion 


formula.) 


bi — pu 
p.. 


Appropriate chordwise integrations of Eq. (32) lead to 
the necessary aerodynamic loads for aeroelastic appli- 
cations. 

The following are some remarks which will be of help 
in practical calculations: 

(1) Inasmuch as the hypersonic small-disturbance 
theory is a steppingstone in the derivation, it is neces- 
sary that r <1. In general, one cannot give an upper 


Py\f y- ) a ‘halal Pee (<2) 21 (x) ) 
= —9 M vyA(x)? My t 
(br) , ( 2 ws is ‘ag ut dx IO + U.. dt T 
y-1. ay dz; . Pe 3 21(0) df(to) } og 
af - | 1 a(s)) nee | (i)! (: U ) Swe a aid 


limit on 7 since this will depend on Mach number and 
the percentage error that one may tolerate. However, 
most pointed wings are sufficiently thin to meet this 
requirement, at least if they are not set at large initial 
angles of attack ap. 

(2) The accuracies of the approximations (22)—(24) 


and (27) are indicated in Figs. 17(a-b) of reference 27. 
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Fic. 3. Low-frequency flutter derivatives for a double-wedge airfoil; comparison of piston theory with shock-expansion method 


However, those figures are not indicative of the errors 
in Eqs. (25)-(26). In general, the latter are slightly 
larger. For better accuracy, the exact transcendental 
shock equations may be utilized with no particular 
difficulty, although they do not satisfy the hypersonic 
similarity rule. On the other hand, if the approximate 
Prandtl-Meyer expansion (27) is replaced by its exact 
counterpart, extensive computational difficulties arise. 

(3) The above two-dimensional derivation may be 
applied with certain reservations on a strip basis to any 
chordwise section of a finite wing. Eq. (32) is inap- 
plicable near a cutoff wing tip; and of course, any strip 
method does not account for aerodynamic span effects. 
However, except for wings with very small aspect 
ratios, very large leading-edge sweeps, spanwise geo- 
metrical discontinuities or rapid modal changes, these 
effects may be safely ignored at hypersonic speeds. 

As a simple illustration of the above, calculations 
were carried out for the usual flutter derivatives,** 
Iy,...L4', My’... , M4’ of a symmetric double-wedge 
airfoil with maximum thickness at 50 per cent chord. 
The reduced frequency was assumed small, so that only 
the leading terms in k were retained. Typical results 
for the moments are shown in Fig. 3, and compared 
with those from second- and third-order piston theory. 

It may be easily verified that the expressions, accord- 
ing to the present analysis or piston theory for the 
quantities RM..M.', k°M..M;’, RM.My', which 
are the ordinates in these plots, obey the hypersonic 
similarity rule—i.e., that they are dependent on Mach 
number and 7 only in the combination JM..r. 

The quantities k°L;’ and k?M;’ (since higher order 
terms in k are neglected in the computations) represent 
steady-state lift and moment coefficients. Another 
important consideration is the location of the aerody- 
namic center, (A.C.), which must be determined accu- 
rately, for flutter speeds are very sensitive to it in prac- 


tical application (cf. reference 31). The A.C. is in 
per cent chord from the leading edge 


A.C. = 100(k?M..M;'/2k?M..L;3’) 


Taking the extended shock-expansion method as a 
standard, it can be seen that second- and third-order 
piston theory start deteriorating at roughly W.r = 0.3 
and 0.7, respectively. Inasmuch as it is simple enough 
to apply, particularly when high-speed computers are 
available, the shock-expansion method is preferred over 
third-order piston theory whenever JV..7 exceeds 0.3. 
An interesting observation is that the quasi-steady 
approximation (utilizing the steady-flow equation, but 
with local effective angles of attack according to the 
unsteady boundary conditions) yields rather close 
agreement for rigid-body motions at moderate values 
of M.r. 


identical. The same cannot be said, however, for 


In the case of pure translation, the results are 


elastic chordwise modes and/or at high values of 
M.t. 

The generalization of the shock-expansion technique 
to steady three-dimensional flows!’ suggests a parallel 
extension of the procedure for airfoils to three-dimen- 
sional bodies. As shown by Eggers, Savin and Syvert- 
son (see also reference 28—30), when J/., > 1 and the 
parameter J/..7 is unity or greater, the flow appears 
locally two-dimensional in planes tangent to the surface 
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Fic. 4. Coordinate systems and notation for a body of revolution 
in unsteady motion. 
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streamlines containing the normals to the body at the 
points of tangency. If A/., times the curvature of a 
given streamline in such a plane is large compared 
with the curvature in a cross-sectional plane, the sur- 
face streamlines may be approximated by geodesics 
originating at the pointed vertex. A more detailed dis- 
cussion is given in reference 13.* The lower limit on 
\/..r depends on body shape; in the above reference, the 
criterion of /..r > 1 has been based on some computa- 
tions on nearly parabolic ogives of revolution with 
maximum diameter at the base. A further require- 
ment for inclined bodies is that the local angle of attack 
a is always much less than unity. For bodies of 
revolution executing small unsteady motions about a 
mean position ay 0, the surface streamlines may be 
assumed to follow the meridian lines. 
Consider the situation of Fig. 4. 
(11), the expression for the equivalent steady-state 


According to Eq. 


body is 
B(x, y, 2) = R(x) — 
- ai(x)f (to + — Vee 
y* 2— 2(X)f | lo = (33) 
satin U./ JS ” 
The boundary condition on the body B = O yields 
(u= U.) 


The pressure on the body is 


= CNG) 
p doI\Pu/ yawes cove \Pw/expansion 


p 


This corresponds to and is governed by the same restrictions as Eq. (28). 
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: v(x, , to) = = + 
—_e dx 
sin y _ f (« + zs ) + 2(x) L SO) | 
dx * U Ue @ \Vaidis/Ve 


(34) 


where v, is the velocity of a fluid particle normal to the 
body and y is the meridian angle. The body-axis 


inclination at the nose is 


é€n(fo) = 
Ox 


dz, l d| ly) 
f(to) + 2,(0) — 
pete l at 


The body described by Eq. (33) is slightly deformed, 
but it is assumed that streamlines still do follow con- 
stant y-lines. By symmetry, this condition is met 
exactly at the windward and leeward sides [W = (7/2) 


and (32/2)|. Following the steps outlined for the air- 


02;(xX, to + XxX) | 


foil, one can set 


py = Px[l + new sin y), 
My = My[l + mey sin y} 


(36a-b) 


where Py and /y are the unyawed values at the nose 
behind the shock, corresponding to @ = @y, and may be 
found, for instance, from conical-flow tables.*# 


2 = 3 7 ° v,(x, W, to) — v,-(0, W, to) 2y/(y¥-1 
a (?") i] + nen sin y} ,! Ae (? - ) ava + men cin »| (x, W 7 ( y } 


Expanding Eq. (37) binomially, retain- 


ing only linear terms in ey, and replacing fo by ¢ — (x, U,) to return to the actual body, one has 


a (P*) 4 (? 
p p \ 


1 - (y¥+1)/(v¥-1) 
) I yA(x)t HE + ? 
hs 


—1\ _ 
= ) My ace | + 


” ] dz, . ; 1 df(t) ' y-l1. as 
yMvy sin y F 1@) + ix) — - + sin WY] ntl + = My A(x) ) + yMxy(mA(x) -— 1) ] X 
dx 2 


l dt 


where 


dR(x) dR(x) ; 
A(x) = : _ (30) 
dx AX \nub 


Lateral loads can be obtained by suitable integrations 
of Eq. (38). 

The pressure expression (38) contains the rate 
parameters 7 and m. Kopal*® has tabulated 7 accord- 
ing to Stone’s first-order solution. It corresponds to 


* Although in principle one may find the geodesics of a body 
of arbitrary cross-sectional shape, experience has indicated that 
difficulties arise near the apex, except for those geodesics in planes 
of symmetry of the body. An application to a body of nearly- 
elliptic cross section revealed that only geodesics with certain 
orientations can emanate from the apex which is a singular point 


(<) f (: x ) ‘ 2(0) df(to) Lo 
— P (o« 
dx Jeo’ U ae? ee 


) 


the quantity designated by 7// in reference 35, while 


m can easily be shown to be 


m l x 1 /n 

M. M.elu, 2 ¥ ) ” 
Ferri*® has criticized this tabulation by pointing out 
that singular points exist around a yawed cone which 
the Stone solution does not account for. He introduces 
the concept of a vortical layer, and finds a more accurate 
result for the flow parameters on the surface of the 
cone and, under the vortical layer, in terms of tabu- 
lated quantities outside the layer. When ey < @y, a 
close examination of Ferri’s correction reveals that the 
vortical layer is very thin and sustains no pressure 


itr 
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difference to O(ey). Eq. (36a) is then valid both inside 
and outside the layer, provided 7 is taken from the 
tables*®® and only lateral loads are being computed. 
The m from Eq. (40), however, gives the sinusoidal 
Mach number variation outside the layer. For the 
corresponding quantity m, at the surface 


Me m 


M M 


[ ce l i é ") (41 
27M v(y — 1)M.3(My/M.)? , p p 


Differences between m, and m increase with M..6y. 

On the basis of the hypersonic similarity rule, the 
tabulations for a 5° cone*!:*® may be reinterpreted as 
functions of M/..@y and may be utilized for calculations 
on bodies of revolutions. In some cases, it is suffi- 
ciently accurate for large M.. and small @y to use the 
following approximate formulas derived from reference 
28: 


ty V1 + [7 + 1/2] (MeO)? X 


AEROSPACE 


SCIENCES—MARCH, 1961 


M.B = V1 + [(y + 1)/2](M.6y)? (42) 
Py/p.o = [1 + ¥(M.0y)*|(My/M,)~?”Y°- (48) 
My _ (2:)\(¥) . | 
uM. \m,/\u (44) 
where, in turn, 
(7 y . 1+ [(y + 1)/2](M.6y)? 
M../ [1 + y(M..6n)*]) 11 + (vy — 19/2] (M.0x)?! 
(45) 
My\? — ' = 
= 4] J 2 
() +( gf Metu \i] * 
M.6y\2 M sy") 
1+1 +. i 
| i Gy 3) = g (6) 


From Savin’s* results for the inclined cone,* one can 
derive after some lengthy algebra 


M 
(2) (M..6n)* Lo 17 
- ; - — > 47) 
\\M, t1 + [(y + 1)/2](M.0n)?} 11 + [(y — 1)/2] (M.6y)?} [1 + v(M.0y)?]h : 
m/M. = —[(y — 1)/2]) V1 + [(y + 1)/2](M.6y)?(Ty/M.)°O (48) 
“—— — 1)/2 4e ; 
m,/M. = —I(y 1)/27] (n/ Mo) (49) provided M.. > 1. Thus, a problem of three-dimen- 
where sional steady or unsteady flow can be reduced to an 
4 unsteady two-dimensional one, in the manner implied 
E +; (A1.0w)" ; | by the simple illustration of Fig. 5. For example, a 
in iL + [(y + 1)/2](M1..6n)? 5 (50) figure of revolution at zero incidence whose shape is 
y¥+5 (M..6n) r = R(x) becomes, for the fluid slab striking its apex 
y¥+1 Vi+ [((y + 1)/2](M..0n)? at ¢ = 0, a cylinder with instantaneous velocity of 
expansion 
The Prandtl-Meyer expansions along the meridian ; 
lines may be done either just outside or inside the v(t) = U at (51) 
dx r=U wt 


vortical layer. For the present case of finding the 
lateral loads only, it is more convenient to do the 
former. 

Although the foregoing discussion has been confined 
to bodies of revolution, the shock-expansion technique 
appears to be extendable to pointed bodies or wing- 
body combinations with rather arbitrary cross-sectional 
shapes.t| For instance, Savin*’ treats the problem of 
a wing mounted on a body of revolution with the 
assumption that the influence of the wings on the 
body is negligible. The unsteady solution can be 
affected along similar lines. Further details are given 
in a fuller report by one of the present authors.*° 


A Variational Approach to Certain Transient- 
Flow Problems 


As pointed out earlier, the piston analogy is valid 
for a slender pointed body at moderate values of M..r 


t See footnote on p. 215. 


When the body is performing small lateral oscillations, 
the same fluid slab will see an expanding and contracting 
cylinder, translating in the z-direction with a velocity 


w(t). The resulting radial velocity at the (now dis- 
placed) piston surface becomes 
. @f ; al 
v,(¥, t) = l + w(t) sin y (02) 
dx r= op t 


To find the airloads on the oscillating body, at any 
given instant of time, it is necessary to consider several 
fluid slabs. Each of these encounters the nose at a 
different initial time fo. 

As an alternate to the difficult method of unsteady 
characteristics,'4 one may obtain an approximate solu- 
tion by a variational method. The appropriate prin- 
ciple is the one given by Bateman,** which reads for 
irrotational motion, and in the present notation 


* In Appendix C of reference 30, the factor Mo’ws in the de- 
nominator of Eq. C.30 should be in the numerator. 
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le 
6/ = ff tple, + |\V¢l2/2] + f(p)} dVdt = 0 
t Vit 


(53) 


where p is density, V is the disturbed volume of fluid 
and f(p) is a function related to the pressure 


bP = pef'(p) — fle) (54) 


g and p are to be varied, with the constraint that d¢ 
must vanish on any part of the boundary of |’ which is 
not impermeable, and the conditions at times ¢; and 
fp are known everywhere in V(t) and V(t). Then J 
is stationary, and the Euler equations yield the usual, 
well-known field equation for ¢. For the range of 
moderate /..7, there is good justification for assuming 
isentropic conditions, particularly for three-dimensional 
pointed bodies. Any desired degree of nonlinearity 
can be introduced into the statement of this problem, 
depending on the extent of the simplifications made 
in the solution of the variational equations. It is not 
difficult to show that, under the condition of isentropy, 
the integrand of Eq. (53) is nothing more than the 
negative of the pressure—i.e., the equivalent form of 


ta 
a i} (—p)dVdt = 0 (55) 
hi V(t) 


It is to be noted that V is finite in the present applica- 


Eq. (53) is 


tions. 

In analyzing piston problems, when the final condi- 
tions are unknown at f, the statement Eq. (55) must 
be modified to allow for the nonzero é¢ at this upper 


limit, *° leading to 


*/o 
6 | f (—p)dVdt — { |pdg|i-. CV = O (56) 
Jt V(t) V (te) 


“hh (y¥ — 1) bh rer? 
F i= — le, + -\W¢l? dVdt + 
J0 J Vit) a.” 2 
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| a Typical Fluid Slob 





=> 





Fic. 5. Body of revolution in axisymmetric supersonic 
flow, seen in still air coordinates as a two-dimensional circular pis- 
ton expanding and contracting 


The variations are to be taken such that: (1) d¢ = 0 
at the outer wave, (2) the boundary condition at the 
surface is unaltered, and (3) at time ¢,, either the dis- 
turbed volume is zero, or the flow is known everywhere 
and the variations are constrained accordingly. Here- 
after, ¢; will be taken to be zero, with no loss in general- 
ity. It is worth noting that the second term of Eq. 
(56) arises because of the lack of a priori knowledge of 
the flow everywhere at time f. This so-called ‘‘tail- 
term” is very similar to the ‘‘“momentum density term” 
of the general Hamiltonian principle when applied to 
structures in transient motions. * 

For isentropic flow, p = /p(p) and both /p, p can be 
expressed in terms of the velocity potential 


v/(y-1 
b/po = [ — [(7 — 1)/a..?| (e, = vel) 
“= 


P/ Po = [ — (vy — 1)/a."] (¢, +<|V¢ :) 


It follows that Eq. (56) reads 





dV =0 (59) 





Eq. (59) is to be solved by assuming a finite series for ¢ which satisfies conditions (1)—(3) above—i.e., by a Ritz ap- 


proach. 


Such a procedure is tantamount to satisfying in some mean fashion (to a degree of accuracy which is depen- 
ying £ 


dent on the number of coefficients in the assumed series) the equation for the velocity potential 


Special forms of the variational principle are 


*le xs(t) q y¥-—1) 
. (y — 1) ae ci 
. | | [ -_-t 2 ( + e?) | xX 
/70 x p(t) Ga” 2 
dxdt + —~- f 
Ga” xp(te) 


for the one-dimensional piston of Fig. (6a), and 


Xs (2) 


{ 


— oa = (Y =— l)eV7¢ Tr [(y —_ 1) ave Ve ? + (O ar), Ve “| + vel (Were | (60) 


(y — 1) ye sa Ni . ; 
| — ms (Gr + Yr° 2) 0g dx = 0 (6la) 
é.* = 


®/o (y = 1) y/(y¥—1) 
6 | f [ _ — (1 + gy?/2 + ¢’ 2 | dydzdt + 
70 S(t) i 


for the two-dimensional piston of Fig. (6b). 


(y — 1) 


he 
.” S(te) a 


6 


~~ y-1) 
(¢ + g,?/2 + ¢? 2) 6 | dydz = 0 (61b) 
t=t2 
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(b) TWO-DIMENSIONAL PISTON 


Fic. 6. One- and two-dimensional pistons in unsteady motion. 


As a first example, consider the simplest case of a 
one-dimensional piston receding at a constant rate 
v, < 0 (fort 2 0). According to the piston analogy, 
this corresponds to a single turning around a corner of 
small angle 6 = v,/U The flow is free of compression 
waves and the outer wave has a velocity equal to the 
speed of sound. The flow field is ‘‘conical’’—.e., the 
velocities and pressures are constant along the charac- 
teristics x/t = constant. 

One series that may then be assumed for the velocity 


potential is 
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Fic. 7a. Pressure on the face of a uniformly receding pis- 
ton vs. piston velocity; comparison of variational and exact solu- 


tions (y = 2). 
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yg = —0,(a.t — X,)f?/2 + 
J cos nzf — | 
A.(Baot — Xp) > A, : - +... (62) 
n=] nT 
where x, = v,¢ and ¢ is the conical variable 
¢ = (a.t — x)/(a.t — X,) (63) 
Several important points should be made here: At the 
outer wave, ¢ = 0, and both ¢ and ¢, vanish irrespective 
of the constants A,. At the piston, ¢ = 1, ¢, = 7,, 


again independently of A,. The first term on the 
right-hand side of Eq. (62) is the nonhomogeneous 
part, which satisfies all the known conditions, whereas 
the second term is the necessary adjustment and is 
termed the homogeneous part. On the other hand, 
the terms containing A,(m = odd) contribute directly 
to ¢, at the piston and, hence, to the piston pressure. 
The A,-terms (~ = even) affect, only indirectly, the 
piston pressure, since they modify the other A,,-terms. 


By taking an increasing number of terms in Eq. (62), 


the convergence of the solution can be established for 
practical applications. For all the examples tried out 
thus far, only few terms were sufficient. The required 
number will depend primarily on the choice of the non- 
homogeneous part and on the distribution functions [in 
this case (cos mmf — 1)| associated with the A,’s. 
Variations of A,, 6A,, which are equivalent to special 
types of variations 6g, satisfy the restrictions (1)—(3) 
stated previously. 

In order to carry out analytically the indicated integra- 
tions, it is necessary to expand the term with the 
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Fic. 7b. Pressure on the face of a uniformly receding piston 


vs. piston velocity; comparison of variational and exact solutions 


(y = 1.4). 
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1 — [(v — 1)/an7] (Gi + 2’ 2)}) ee oy 
1 — (1/a.*)(¢: + ¢2?/2) + [(2 — y)/2] X 
(1/a..*)(¢: + ¢27/2)? +... (64) 


When the approximation Eq. (64) is used in the vari- 
ational equation, there results a set of simultaneous 
nonlinear algebraic equations which may be solved, for 
instance, by iteration to yield the unknown coefficients 
A,. Once these coefficients are determined, the piston 
pressure follows from Eq. (57): 


l vp 
P itie~ ee a 4+ 
p 2a 
2 Up | wie 
(: - *) } 2 A,$ (65) 
T a n=1,3,5 nN { 


The pressure on the face of a uniformly receding 
piston versus a dimensionless piston velocity is shown 
in Figs. 7a-b. In the first case, y is taken purposely 
to be 2, thus avoiding the truncation error in Eq. (64), 
and concentrating on the convergence of the procedure 
by taking NV = 1, 3 successively. Also shown on the 
same plot is the known exact result for y = 2. Even 
the one-term approximation yields good agreement. 
The example is repeated for y = 1.4, retaining two and 
three terms of the binomial expansion, in Fig. 7b. As 
expected, the inclusion of the third term results in 
appreciable changes at the high values of vp, roughly 
above |v,/a.| > 0.5. 
of v, is due to the one-term approximation, similar to 


The discrepancy at lower values 


the situation in Fig. 7a. 

If the receding piston has a velocity varying with 
time, the flow would no longer be conical, and a velocity 
potential of the following form may be assumed: 

l 
y — 5 va t — X»)F* + a.(Gut — Xp) X 


I | Vu f m—1 
( c3 ) 2, & ( ) +...]| (66) 
3 t te 


Note that the ¢-distribution functions are here taken to 
, which satisfy all the required 


be (¢7/2 — ¢°/3), 
conditions, rather than the functions (1/27)(cos nf — 
1) of Eq. (62). Of course, the latter could have been 
used equally well. 

The ideas set forth above can be carried over to 
treat two-dimensional piston problems. Furthermore, 
the piston motions need not be of the receding type. 
In dealing with two-dimensional expanding pistons 
having v, > 0, two modifications may be necessary if 
v,p/a. is not very small. This is because the shock 
velocity v, may be substantially greater than the speed 
of sound, and the particle velocity just behind the shock 
is not zero as in the case of that behind an expansion 
wave. 

To illustrate, consider the piston problem analogous 


to a circular cone at zero angle of attack. A suitable 
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Fic. 8. Pressure coefficient on circular cones at zero angle of 


attack vs. 7..06v; comparison of variational solution with those 
according to references 16 and 34. 


form for ¢, expressed in terms of the polar coordinate 
r, which satisfies all the requirements of the variational 
principle, is 
gy = —v,(vst — R,p)f?/2 + aa(vst — Rp) X 
— 2 /6 +3 /2\91 9 ' - 
tAi(S — £7/2) + Ao(g?/2 — £2/3),4+5...5 (67) 
where 
¢ = (vt — r)/(vst — Rp) (68) 
At the outer wave ¢ = 0, ¢ = —Aya, + 0. The 
variational principle yields N equations, one for each 
A The additional unknown 2,, 


n* 


stant during the variation, is obtainable from the 


which is kept con- 


Rankine-Hugoniot jump conditions across a weak 


~y l 
[1 — a A, fe (69) 


In Fig. 8, the results of such a calculation are pre- 


normal shock 


sented and compared with the exact solution of Kopal** 
for a 5°-cone and with that of Van Dyke.'® It appears 
that the assumption v, = a.. is tolerable up to J/..@y 
0.5. On the other hand, the A,-term must be retained 
even for the low range of \/..6,. This allows for the 
particle velocity just behind the shock to assume a non- 
zero value. 

A similar calculation for a parabolic ogive with three 
terms of a suitable series was carried out at a low value 
of M.6y. 
slender-body theory in Fig. 9. 

The foregoing examples have been confined to axi- 


The results are compared with linearized 


symmetric flows. Two simple situations of nonaxi- 
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symmetric flow have also been analyzed. The first 
involves an elliptic cone at zero angle of attack, while 
the second considers a circular cone at a small incidence. 
To facilitate the evaluation of the involved integrals, 
certain coordinate transformations were utilized. The 
results for the elliptical cone are presented in Fig. 10 
and compared with Van Dyke’s second-order solution.!! 
The agreement is deemed quite satisfactory in view 
of the fact that only two terms of the assumed series 
were taken. As for the 5° semi-apex circular cone at 
M., = 4.86, the pressure distribution due to a small 
angle of attack turns out to be 


(Ap). = —K(p. asin p) 


with A = 5.596. This result compares reasonably 
with the exact value of 5.435 from cone tables.*® The 
corresponding coefficient from slender-body theory is 
K = 5.786. 

Although the foregoing examples have dealt with 
steady-flow problems, the variational approach need 
not be confined to such cases. The extension to un- 
steady motion is discussed in detail in reference 40. 
Briefly, the steps are quite similar to those involved 
in the case of a slightly yawed body in steady flow. 
If the unsteady velocities are small, the unsteady part 
can be linearized, and the analysis is then a perturba- 
tion solution for the time-dependent motion over the 
main flow field. Moreover, results for various types 
of motion may be superimposed. 

The success and the practicability of the variational 
approach depend on certain factors, the most important 
of which may be summarized as follows: 

(1) Judicious choice for the form of the assumed 
potential. It is desirable to have the nonhomogeneous 
parts as close to the actual solution as possible. To this 
end, linearized or other approximate solutions, if avail- 
able, can serve as valuable guides. 
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Fic. 9. Pressure distribution along a parabolic ogive of rev- 
olution at zero incidence; comparison of variational and linearized 
slender-body solutions. 
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Fic. 10. Pressure distribution along periphery of an elliptic 
cone at zero incidence; comparison of variational! solution with 
that from second-order theory of reference 11. 

(2) Accurate integration techniques for the double 
and triple integrals of the variational equation. For 
almost all of the examples above, the integrations were 
carried out exactly. For more complicated situations, 
it may become necessary to resort to numerical pro- 
cedures. Experience has indicated that extreme cau- 
tion must be exercised in these integrations, since the 
matrices associated with the simultaneous equations 
for the unknown coefficients are not always well- 
conditioned. 

(3) The calculations are quite lengthy and, except 
for the simplest problems, undoubtedly require high- 
speed digital computers to render the method practical. 
With this in mind, the procedure must be highly sys- 
tematized, an objective which will require some further 


research. 


Application to a Static Aeroelastic Prob!em 


The two methods for predicting airloads described 
in this paper are now applied to static divergence of a 
cantilevered ogive of revolution. It is recognized at 
the outset that divergence of such configurations is not 
necessarily important from a present-day design stand- 
point. Nevertheless, it serves as a good illustrative 
example. 

Consider this ogive to be of length / and radius 


oe x If/x\*| " 
ne off 1). ao 


where 6y is the semiapex angle and x is the streamwise 
coordinate with origin at the apex. The body is canti- 
levered at its base, and it is assumed to behave struc- 
turally as a simple beam with bending stiffness 
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El (x) = (EI)o(x/l)’ (71) 


Here (EJ) is the bending stiffness at the root. This 
problem is the exact counterpart of that for a very low- 
aspect-ratio delta wing treated by Martin and 
Watkins."! 

The aerodynamic loads are computed for \/.0y => 
0.8 by the shock expansion method, and for A/.@y = 
0.5 by the variational procedure. Originally, the 
iterative solution proposed in reference +1 was adopted. 
It was discovered that the resultant slopes of the de- 
flection z:(x) of the body axis were rather insensitive 
to changes in \/.@y and that these could be approxi- 


mated satisfactorily by 
dz;(x)/dx = (constant)[1 — (x//)?| (72) 


Eq. (72) was subsequently used for all the calculations. 
The stability boundaries thus computed are shown in 
Fig. 11. The ordinate represents the dimensionless 
stiffness level below which the structure becomes stat- 
ically unstable. The boundary obtained from slender- 
body theory is given for comparison. The curve from 
the shock-expansion method starts to dip at \/..0y ~ 1; 
the result below this value is not expected to be reli- 
able. The single point at 1/.@, = 0.5, obtained from 
the variational method, fairs in well with the shock- 
expansion curve. As_ slender-body theory should 
apply at .\/..0y near zero, the curve is completed some- 
what empirically by passing it through this point. 

To get some idea of the stiffness needed to avoid 
divergence under a given set of flight conditions, con- 


sider a body with 6y = 0.2 radius at 1/.. = 10. One 
finds 
(ET) “ 
— & 0.86 (73) 
q Oy7/3 
With typical values of / = 30 ft., E = 10’ psi, and 
dynamic pressure = 2,000 psf, the required moment 


of inertia at the root is J) = 800 in.4 Assuming this 
due entirely to an unbuckled monocoque skin of uni- 
form thickness, it corresponds to a gage of 0.005 in. 
Since entry vehicles, for example, generally rely on 
internal structure supplying stiffnesses equivalent to 
much thicker skins, one would not expect much trouble 
with divergence even at the rather high peak dynamic 
pressure assumed here. Further applications must be 
made to assess the influence of aerothermal effects, and 
to determine whether a similar conclusion holds for 
vehicles being boosted out of the atmosphere. 


Concluding Remarks 


The principal objective of this paper is to demonstrate 
two approximate aerodynamic theories, suitable for 
use in the analysis of steady and unsteady aeroelastic 
phenomena on a wide class of pointed and slightly 
blunted configurations, at values of M.7 typical of 
atmospheric entry. Both techniques are believed 
adaptable to high-speed digital computers and, there- 
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fore, suitable for synthesis with larger programs for 
calculating dynamic loads, stability, flutter eigenvalues 
and the like. Of the two, the variational procedure 
is the more complicated and requires the greater de- 
velopment effort. Its range of applicability, as pre- 
sented here, is more limited than that of the shock- 
expansion method. Nevertheless, it is considered to 
hold great promise of extension to other situations of 
practical interest, especially if the restriction to isen- 
tropic flow of a perfect gas can be removed. 

Since no direct or indirect experimental evidence can 
yet be offered to establish the accuracy of airload or 
aeroelastic computations, more exact theory furnishes 
the only source of comparison. Additional examples 
of such checks and more detail on both methods will be 
found in the extended report.” A pressing need exists 
for more precise and comprehensive measurements 
of airforce distributions at genuinely hypersonic Mach 
numbers. This remark applies particularly to unsteady 
loadings at appreciable reduced frequencies, in view of 
the desirability of ascertaining, among other things, 
how widely the quasi-steady approximation can be 
employed. 

The example of body divergence given in the fore- 
going paragraphs was intended more to illustrate the 
feasibility of the new theories than to identify an aero- 
elastic problem of pressing importance to the designer. 
Nevertheless, it should be noted that the unbuckled 
cylinder is a very efficient bending structure, and that 
flatter wing-like vehicles may be susceptible to this 
sort of chordwise instability (as already found by Mar- 
tin and Watkins*'). Hypersonic and entry aircraft 
will probably pass through a period of configuration 
refinement, during which a broad variety of shapes 
and sizes will be examined and tested. It is, therefore, 
premature to generalize about the relative seriousness 
of stiffness and aeroelastic requirements in this flight 
regime. Only experience and time will determine 
the true utility of aerodynamic tools such as those 


reported herein. 
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The Buckling of Thin Spherical Shells 


F. J. MURRAY* anp FRANK W. WRIGHT** 


Duke University 


Summary 


Precise numerical solutions are obtained for the energy and 
equilibrium equations given by Th. von Karman and H. S. Tsien 


for thin spherical shells under uniform external pressure. The 
method of solution employs a power series expansion to move 


away from the singularity at the initial boundary, and a step- 
by-step integration procedure throughout the desired angle open 
ing of the shell. The ‘‘small angle’’ assumption is eliminated 
without introducing trigonometric computation. For fixed ra- 
dius and thickness, the various equilibrium configurations of the 
shell are obtained as a function of two parameters, one corre 
sponding to external pressure and the other corresponding to the 
curvature of the shell at the central point 

The solutions indicate clearly the behavior of the shell in both 
the pre-buckled and post-buckled states, and yield distinct 
values for the upper and lower buckling loads. These critical 
loads are determined on the basis of the load-deflection char- 
acteristics. The results obtained depend critically on terms 
which are neglected in the classical theory by the process of 
linearization 

The same computing procedures were applied to the mathe- 
matical formulation of Keller, Greenberg, and Reiss and yield 
an improved resolution of the phenomena associated with 
buckling 

The calculations were performed on the IBM 650 Electronic 


Data Processing Machine 


Symbols 
a,b,c,d = quantities defined by Eq. (9) 
B = quantity related to the angle opening of the 
shell—see Eq. (13) 
Cy Es = power series coefficients 
E = Young’s modulus 
, = derivative function 
h = step interval in the numerical integration pro- 
cedure 
H, & = dimensionless energy quantities 
= curvature of the shell 
p = external pressure 
P = dimensionless pressure term 
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q = parameter corresponding to the curvature of the 
shell at the central point—-see Eq. (30) 
r = horizontal distance from the axis of symmetry of 


the shell 


R = radius of the shell 
; = parameter corresponding to the stress at the 
central point of the shell-—-see Eq. (59 
t = thickness of the shell 
u = dimensionless pressure term—see Eq. (1 
U3, Uo = upper and lower buckling loads respectively 
V,s, Z = auxiliary variables in the numerical integration 


procedure 
= radial deflection of the shell 


W = total energy of the shell 

v, y = transformation variables—see Eq. (8) 
a = axis of symmetry of the shell 

a = polar angle 

B = semi-angle opening of the shell 

7 = stress function—see Eq. (54 

6 = deflection of the shell at the center 
7] = angle of inclination of a meridian 
0 = slope function—see Eq. (60) 

r = geometric parameter 

A = semi-angle opening of the shell 

v = Poisson’s ratio 

I] = total energy of the shell 

p = geometric parameter 

o = stress quantity—see Eq. (54 

¢ = error function—see Eq. (52 


Introduction 


A FUNDAMENTAL DISCUSSION of the theory of the 
buckling of thin spherical shells is given in a paper 
by von Karman and Tsien! which derives an expression 
for the energy and the equation of equilibrium for the 
shell when subjected to uniform external pressure. 
Th. von Karman and Tsien! obtain an approximate 
solution for the minimum buckling load by the Ray- 
leigh-Ritz method. In the present paper, techniques 
are prescribed which yield precise numerical solutions to 
the equation of equilibrium, and these results clearly 
describe the buckling phenomenon known as _ the 
‘“Durchschlag”’ or “‘oil can’’ effect. The major ob- 
jectives of this investigation are to permit a comparison 
between the precise results and those given by the von 
Karman-Tsien! approximate method, and to acquire 
an appreciation of the nonlinear character of the prob- 
lem. A comparison of this numerical procedure is also 
made with the method of finite differences given by 
H. B. Keller and E. L. Reiss.’ This comparison is 
made using precisely the same mathematical formula- 
tion as was first given in a paper by Keller, Green- 
berg, and Reiss.® 

The equation of equilibrium given by von Karman 
and Tsien! is a nonlinear second-order differential equa- 
tion having a singularity at the origin. The method of 








3 Axis 











Fic. 1. Geometry of the spherical shell. (Note that # is nega- 
tive. Fig. 3 of reference 1 shows @ as a positive angle and is in- 
consistent with the formulas of that paper. ) 


solution uses a step-by-step numerical integration on 
an automatic sequence calculator. The singularity in 
the equation is handled by a power series expansion, 
but this power series is valid only over a relatively small 
fraction of the region of interest. In specific cases, for 
example, this fraction is only one seventh. Thus, the 
power series expansion is used simply to carry the solu- 
tion away from the singular point, and thereafter the 
step-by-step integration procedure is used. An alterna- 
tive procedure using the Runge-Kutta method for start- 
ing is also possible and was utilized in the computation 
of the Keller-Greenberg-Reiss® equations. 

The immediate output of the computation is a two- 
parameter family of solutions of the equilibrium equa- 
tion for a shell with a given radius and wall thickness. 
Specific situations are then described by imposing 
boundary conditions which yield a one-parameter 
family of solutions. A one-parameter set of solutions 
of this type describes the possible configurations of the 
shell for various pressures. 

It was pointed out by Friedrichs? that the assump- 
tion of vertical displacements made by von Karman 
and Tsien! is restrictive, and a revised system of equa- 
tions was published by Tsien.* In the present dis- 
cussion, numerical solutions are obtained for both the 
von Karman-Tsien equations and the Keller-Green- 
berg-Reiss® formulation of the Friedrichs correction. 
However, the latter formulation involves a number of 
simplifications, including the small-angle assumption. 

As mentioned above, the Rayleigh-Ritz solution is 
given by von Karman and Tsien.' In addition, there 
have been papers based on the classical linearized 
theory, brief accounts of which are presented by 
Timoshenko‘ and Kaplan and Fung.’ Recently, there 
have been investigations based on the use of a power 
series expansion, examples of which are those of Reiss 

sreenberg, and Keller,® and Weinitschke.’ 

The Rayleigh-Ritz method involves a simplification 
of the equilibrium equation by assuming small angles 


and dropping terms. The present investigation shows 
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that the Rayleigh-Ritz solution is not a solution of the 
equation even in an approximate sense. It is also 
true that the solution of the complete differential equa- 
tions depends critically on terms which are omitted in 
the classical theory. These terms are either negligible 
or they introduce singularities into the solution which 
are ignored in the classical equations. The power 
series investigations do not yield the full picture of 
buckling phenomena as described above. In _ these 
approaches the critical buckling load is associated with 
the divergence of the series, and the extent to which 
these two phenomena are related is not obvious to the 
present authors. A complex singularity in the solu- 
tion will limit convergence, as will real singularities. 
In a more recent paper Weinitschke® modifies the power 
series method by applying analytic expansions to extend 
the radius of convergence. Keller and Reiss’ abandon 
the power series method altogether in favor of the 
method of finite differences. The solutions given in 
references 5, 6, 8, and 9 are restricted by convergence 
difficulties to a limited range of shallow shells, and 
the results yield values for the upper buckling load 
only. The solutions of reference 7 yield values for the 
upper and lower buckling loads but are restricted by 
certain simplifying assumptions such as the use of small 
angles. 

The basis for the present numerical procedure is the 
well-known step-by-step method of integrating differ- 
ential equations. As mentioned above, some modifi- 
cation is essential in order to handle the initial singu- 
larity which appears in both systems treated in the 
present paper. There are two basic advantages in the 
step-by-step procedure which are obtainable simply at 
the cost of additional computation. The first is that 
terms need not be dropped. The second is that any 
accuracy is ultimately attainable. 
yields, in the Keller-Greenberg-Reiss® equations for 
instance, a finer resolution of the buckling phenomena, 
which would otherwise be obscure (see the last section 
of this paper). On the other hand, the dropping of 
terms in the Keller-Greenberg-Reiss® derivation may 
have critically affected the phenomena by eliminating 
singularities such as are present in the von Karman- 
It seems very desirable to apply our 


The higher accuracy 


Tsien! equations. 
numerical method to a complete, unsimplified system 
of equations. These procedures also incorporate a 
method of eliminating the small-angle assumption 
without introducing the need for trigonometric compu- 
tation. 

The parameters which appear in the computation are 
u, the normalized pressure term, 8, one half the solid 
angle subtended by the spherical segment (see Fig. 1) 
and q, a quantity corresponding to the curvature of the 
shell at the apex. The parameter uw is defined by the 
equation 


uw = pR/2Et (1) 


where p is the external pressure, XR is the radius of the 
spherical segment, ¢ is the thickness, and E is Young's 


modulus. The computation yields a fundamental 
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relation between 1, g, and 8, and all results are derived 
from this relation. It would be highly desirable to 
verify this relation experimentally since it is essentially 
independent of initial imperfections in the test speci- 
mens which may lead to discrepancies in the buckling 
pressures. The angle 8 and the ratio K/t would be 
prescribed by the given shell, « would be the inde- 
pendent variable, and g the dependent variable. The 
curvature g should be measurable by optical means if 
the shell is polished to a mirror finish. 

The numerical computations have been carried out 
in complete detail for one case, that of a clamped 
spherical shell with a ratio R/t of 362.5. The results 
of the von Karman-Tsien case yield values correspond- 
ing to the upper buckling load m, and the minimum 
buckling load w, as follows: 


u, = 0.001126 wo = 0.001077 


The critical values of u given above are independent of 
the included angle 28 when it is in the approximate 
range 0.388 < 28 < 7a radians, i.e., when the buckled 
region is confined to a solid angle of approximately 
0.38 radians. The standard eight-decimal floating 
point facility available in the machine was utilized 
and even this introduced limitations on the precise 
determination of the size of the dimple and also of the 
deflection of the shell at the center. The effect of this 
limitation and the character of all results obtained are 
discussed below. 

The step-by-step integration procedures were also 
applied to the equations given by Keller, Greenberg, 
and Reiss.6 The numerical results agreed with those of 
Keller and Reiss,’ but considerably more phenomena 
were observed and are described below. However, the 
mathematical simplifications used in the derivation of 
the Keller-Greenberg-Reiss® equations indicate a need 
for caution in regard to the significance of these phe- 
nomena. Some of the terms dropped, for example, may 
have produced singularities which would nullify the ob- 
servations. This again indicates the desirability of a 
more complete discussion. 


General Discussion of the von Karman-Tsien 
Energy Equations 


Given the geometrical configuration of a spherical 
shell as shown in Fig. 1, it is desired to obtain the load— 
deflection characteristics and to determine the minimum 
external pressure which will cause buckling. The von 
Karman-Tsien! system of equations is based on energy 
considerations which yield equilibrium configurations 
of the shell for variations in pressure. The assump- 
tions made in reference | are as follows: 

(1) The deflection is rotationally symmetric. 

(2) The deflection is vertical (i.e., parallel to the axis 
of symmetry). 

(3) The effect of lateral contraction is neglected 
(1.e., Poisson’s ratio is assumed to be zero). 

(4) The solid angle of the shell is small. 

Assumptions (1) through (3) are retained in this 
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discussion. The small-angle assumption, however, 
is avoided. 
The equations of interest correspond to Eqs. (3), (4), 


and (5) of reference 1 and are summarized as follows: 


W, E t [ COS @ 
— = Jo — 1) sinada (2) 
Rx RJ, \cos 6 ; 


Ws l YS "31 /cos 0 . 
= E } 6’ -- 1 + 
R*z 1? Js COS a 
sin 6 F.. 
: — | snada (3) 
sin a 


W; °8 
Rn p | sin? a(tan @ — tan a) cos ada (4) 
ie ‘ 


where /V’, is the strain energy due to the extension of 
the shell elements; ]V, is the strain energy due to bend- 
ing; and |W, is the work done by the pressure in de- 
forming the shell. 
The total energy of the system, WW, is the sum of 
Eqs. (2), (3), and (4) multiplied by the term R*r. 
The deflection at the center is given by Eq. (12) of 


reference |: 


8 
i { (tan a — tan @) cos ada (5) 
e 0 


The equation of equilibrium as given by Eq. (7) of 


that reference is: 


_f sin a cos a cos a 
2k: tan 6 — 1] 


“~R cos 8 cos 6 / 
l k t\3 | (= 4 
5 cos 6\ tan- - 
6 RJ tL sin a T _ 
cos” 6 de 
(2 tan? a + 1) 

cos @ da 

sin 6 cos 6 tan a (7) cos? 6 tan a d’6} 

COS @ da COS a dat | 


p sin® a cos a sec? 6 = 0 (6) 
The boundary conditions are 
@=Oata=0; 6=Bata=8 (7) 


Without the aid of automatic computation, this 
equation is virtually impossible to solve. As an alter- 
native, von Karman and Tsien assumed a special form 
of displacements and Rayleigh-Ritz 
The equations were further simplified by 
these 


applied the 
method. 
assuming small angles. As a consequence of 
simplifications, the nonlinear terms were dropped. 
The results of the present paper show, however, that 
such terms cannot be neglected and that neither of the 
above simplifications yields solutions to Eq. (6). 


Method of Solution 


The small-angle restriction and the use of the Ray- 
leigh-Ritz method can be avoided by a direct computa- 
tional attack on the equivalent of Eq. (6). The 
method of solution described in this paper uses a power 
series to move away from the singularity at the central 
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point and thereafter employs a step-by-step integration 
procedure throughout the desired range of 8. An 
attempt was made to solve this equation solely by the 
method of infinite series similar to that described 
by Reiss, Greenberg, and Keller® and Weinitschke,’ 
but serious convergence difficulties appeared. This 
method was therefore abandoned in favor of the one 
described herein. 

It is desirable for computational reasons to eliminate 
the trigonometric functions in Eq. (6). This is ac- 
complished by making the following change in vari- 


ables: 
x = tan (a/2) y = tan (6/2) (8) 
Now let 
a=1-22 cm1-y 
: - > (9) 
b=1tx% d=14y% 


Therefore the derivative expression corresponding to 
the total energy given by Eqs. (2), (3), and (4) becomes 


_ Ww ad \?x 
” “ja, * 5 ~ ) et 
1 ft} (b*c ' by t 
3 () | (o "™ ) ” (; x ) 1b? i 
32u (‘ y— :)  % 10) 
c bt 


where the prime (’) indicates differentiation with 
respect to x, and Eq. (5) yields 


6’/t = 4(R/t)[x — (a@/c)y|(1/b?) (11) 


Eq. (6) may be written in the form 
P y ad 4 la) 1“ (7) a*d' (= ) 
Y ~ x Be xb~ ~ Nt) bet L? bc 7 i) + 
d a? t+a 
ux. | + 4x — — 2x y’ + 
b ch} ab : 


(2£4) on 
2y (y)* (12) 
cd 


The boundary conditions corresponding to Eq. (7) 
are 
y=Oatx=0; y= Batx=B 
where B = tan (8/2) (13) 


Inspection of Eq. (12) shows that for 0 < x, y < 1, 
the terms on the right-hand side are well behaved, 
while the left-hand side contains the singular point at 
x = 0. The singularity may be handled by manipu- 
lating the left-hand terms as follows: 


ie yad{ 1a _le (14) 
, xbelxb yd-~ 


The bracketed term can be expressed as 
(d/dx)[F(x) — F(y)] (15) 

where F(x) = log [x/(1 + x?)]t 
F(y) = log [y/(1 + y*)]f 


(16) 


Hence, Eq. (15) becomes, upon substituting Eqs. 
(16), 


yb 
(d/dx)[F(x) — F(y)|] = —(d/dx) tog (: | (17) 
x « 


and Eq. (14) can be written as 
, , vad d , (2 i) " 
yt x bc dx L 8 \xd of tie 
ad? d ? ') ad d (?) yad (; 
cb? dx \xd] ~~? bc dx \x x cdx \b 


Thus the equivalent of Eq. (6) in terms of the quanti- 
ties x and y is obtained: 


ad d (”) (=) ee : (= ) 4 
" pa = 48 , ee 
+ bc dx \x ; t b'c4 | \be 
d d? “ (° aa *) : 
rae ea ee 
2+¢ yad|y x 
2y (y’)? + 2° — yy’ — 
° ( cd ) wy + x be ; , | ie 


This is suitable for manipulation with infinite series. 
The quantity y/x is undetermined at the origin and its 
limit as x approaches zero, q, is considered as a param- 
eter throughout the remainder of the computation. It 
has the character of a constant of integration in our 
discussion. Geometrically, it corresponds to the curva- 
ture of the shell at the central point, as explained at the 
end of this section. 

Eq. (18) is reduced to a system of two first-order 


~ 


equations by introducing auxiliary variables. The 
first equation can be obtained as follows: let 


2 = y’ + (y/x)(ad/bc) (19) 


and v= xy (20) 


Then, by differentiating Eq. (19) and substituting in 
Eq. (18), we obtain 


R\? a*d' {4d 2y aN 
2’ = 48 yun (x? — y*)¢ + 


t bec* b bc 
nad jad 2 d ) | \ P| 4 
ab \be Lb (2 + @) <i (4 + a)y ( 


(2 +)c) vad 
2Qyy’ yo + (21) 
oi cd x be 


The second equation may be obtained by differenti- 
ating Eq. (20) and substituting the expression for y’ 
from Eq. (19). The relations of Eq. (9) are also used. 
Hence, 


v’ = ax + (2y/bc)(x? — y?) (22) 


y’ = 
Thus, Eqs. (21) and (22) constitute the pair of first- 
order differential equations required for the solution of 
Eq. (12). 

Eq. (21) expresses 2’ as a function of x, y, and y’. 
However, y can be expressed in terms of v and x and, 
by Eq. (19), y’ can be expressed in terms of x, 2, and 2. 
Thus, effectively, Eq. (21) can be considered to be in 


the form 
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z’ = f,(x, v, 2) (23a) 


v’ = fo(x, v, 2) (23b) 


It should be pointed out that the possibility of a series 
solution is apparent directly from Eq. (12) if the left- 
hand terms are written: 


y” = (a/b) [(xy’be — yad)/x?] (24) 


However, the computation is made considerably easier 
by the above manipulation. 

That qg corresponds to the curvature of the shell at 
the origin can be seen as follows: 

We introduce Cartesian coordinates (7, 3) where 3 
is the axis of symmetry and r is perpendicular to the 4 
axis. The shell is obtained by rotating a curve 4 = 


f(r) about the 4 axis (see Fig. 1). Then r = R sin a, 
and by Eq. (8) 
x = (r/2R)(1 + x?) (25) 
The slope of the curve 4 = /(r) is given by 
da —2 tan (0/2) —2y 
= —tand@ = = = > (26) 
dr 1 — tan? (0/2) l— y’ 
and hence, 
y (da /dr)(1 — y?) = 
—-~ = —R a (27) 
x r(1 + x°) 


Due to the assumed symmetry of deflection, 


da/dr|\,.0 = O 


and since lim y = 0, it is clear that 


x—0 


lim (y/x) = —R(d*3/dr*) (28) 


x—>0 
The curvature K of the curve 4 = /(r) is given by 


d*3 dr? 


K = 7 (29) 
[1 + (da/dr)?|” 
Atr = 0, da/dr = 0 and therefore we may write 
g = lim (y/x) = —R(d*%3/dr?) = —RK (30) 


x—0 


Thus, if K < 0, then g > 0 and the central portion of 
the shell is concave downward and is representative of 
the pre-buckled state. If A > 0, then g < 0 and the 
central portion of the shell is inverted as in the post- 


buckled state (see Fig. 1). 


The Infinite Series 


The method of infinite series provides a means of 
moving away from the singularity at the origin and also 
of obtaining the starting values for the step-by-step 
integration process. 

Since the deflection is assumed to be symmetrical, 
6(a), or its equivalent y(x), is an odd function (see page 


19 of reference 1). Thus, the series for y takes the form 


y= Cx t+ Co? + Co + Cr" +... (31) 
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Eq. (31) satisfies the first boundary condition 
(y = 0 at x = 0) while the derivative of y (Eq. 31) 
evaluated at x = 0 is C, and is undetermined. 

Dividing Eq. (31) by x, the quantity y/x evaluated 


at x = 0 is also equal to C, which is undetermined. 
The coefficient C, is the parameter g previously defined 
Since v = xy (Eq. 20), the series for v is obtained from 


Eq. (31) and is given by 


| | 8 


v= Cx? + Coxt + Csx6 + Cr Yt : (3: 


to 


and v’ is obtained by differentiating Eq. (32). 
The quantity z may be obtained by expanding Eq 
(19), which yields 
z= Ko + Kox? + Kyx' + Kort +... (33) 
and z’ is obtained by differentiation of Eq. (33). The 
coefficients C and K are obtained by successive differen- 
tiation of Eqs. (18) and (19), respectively, and by 
evaluation of the results at x = 0. These coefficients 
are polynomials in g, and in this solution the first four 
nonzero terms of the series are considered. These 


are shown in Appendix A. 


The Step-by-Step Integration 


A brief description of the method of integration used 
has been published by F. J. Murray.'° This method is 
a modification of the Milne'' method and requires 
knowledge of four previously computed values of the 
derivatives taken at equally spaced values of x. 

Consider the system of differential equations given 
by Egs. (23): 

z’ = f(x, v, 2) (23a) 
v’ = fo(x, v, 2) (23b) 
Let x, = nh, where m is an integer. The method of 
integration assumes that the various quantities v, v’ 
z, and z’ are known at Xy_4, Xy»—3, Xn—2, X»—1, and obtains 
To initiate the process, the power series is 
X_2, X11, and 


them at x,,. 
used to determine these quantities at x_,, 
xo. The values for //’ and 6’/t at the corresponding 
points are given by Eqs. (10) and (11) respectively. 
Once these values are determined, the solution is con- 
tinued by successive open steps throughout the range 
of x. At prescribed intervals, a closed step is taken 
to check the accuracy of the previous open step. The 
closed-step values are not used to correct the previous 
open steps but a comparison of the two is used to deter- 
mine the appropriate interval of integration. As dis- 
cussed by F. J. Murray,'® the accuracy of the complete 
solution can be verified by using a smaller interval of 


integration. 
The equations for the open step are 
. Of + 37f 
e (0) — o i(—Di}.e—-4 TT O4f]1.s —_ 
Zn Z.1 + 24 ; 
IM ne + 5dfi.n-1) 
(34) 
h 
~~ ee a (—Ofe ,-4 + 37fe., — 
n Un—! 24 
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Fic. 2. Typical load—deflection characteristics. 
where /: is the interval of integration and the second 
subscript of /1,,-: implies that the value of /; is taken 
fOr Xp_1, Zn_1, ANd V,_1. 

The corresponding formulas for the closed step are 


h 


Zn” = Zn-1 a 24 (fi.n ze fin 2 + 


19/1 2-1 + fi.) 


5 a 
Un” = Vy} + 24 (fe,n —" OSe.n 2 + 


19fo.n-1 + Of2.n) 


Eqs. (10) and (11) may also be expressed as functions 
of x, 2, and v by use of the relation y = v/x and Eq. 
(19), i.e., 

H’ = f;(x, v, 2) 6’/t = fi,(x, 9, 2) (36) 

We can then integrate /7’ and 6’/t by formulas analo- 

gous to Eq. (34) in a step-by-step manner and obtain 
the values of these quantities at each x, 

The quantities y and y’ are immediately expressible 


in terms of x, z, and v. 


The Technique of Solution 


The system of equations discussed previously yields 
three basic solutions which can be expressed by the 


following relations: 


y = V(x, g, u, R/t) (37) 
6/t = D(x, gq, u, R/t) (38) 
H = H,(x, q, u, R/t) (39) 


Relation (37) represents the solution of the second- 
order differential equation of equilibrium given by Eq. 
(12). Relations (38) and (39) correspond to the solu- 
tions of Eq. (10), the ratio of the total deflection at the 
center to the thickness and Eq. (11), the total energy 
of the system, respectively. For a given R/t, u is 
fixed in the computation and q is varied through a given 
range. This procedure is repeated for variations of 1, 
which yields a two-parameter family of solutions for 
the desired ranges of g and 1. 


In this procedure, the solutions of interest are those 
which satisfy the second boundary condition, i.e., y = B 
atx = 8B. This can be expressed by the relation 


B — Y(B, q, u, R/t) = 0 (40) 


which is obtained by substituting B for x and y in Eq. 
(37). The quantity B is related to the solid angle of the 
shell by Eq. (13). The solutions, Eq. (40), yield the 
relation 


B = F(q, u, R/t) (41) 


which can be expressed graphically as a family of curves 
with one graph for each value of R/t, with various curves 
corresponding to the variations of . 

It is desired to determine the relations for 6/t and 
IT as functions of B, u, and R/t. This is done by elimi- 
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after buckling. 
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nating g from Eqs. (38) and (39) in the following 
manner. 
From the graphs, Eq. (41), we can form the relation 


q = Q(B, u, R/t) (42) 


Substituting Eq. (42) into Eqs. (38) and (39) the 
following expressions for the deflection and energy 
are obtained: 


é/t = D,[B, O(B, u, R/d), u, R/t] = 
D.(B, u, R/t) (48) 


H = H,{B, Q(B, u, R/t), u, R/t] = H(B, u, R/t) (44) 


These relations may then be used to determine the 
deflection and energy as functions of « for any spherical 
segment given B and R/t. These relations can be ex- 


pressed as follows: 
6 ‘|p -—@ = Ds(u) (45) 
R/t = R/t 


i) _— = Au) (46) 
R/t = R/t 

Relation (45) permits a graphical representation of 
the load—deflection characteristics for a given spherical 
shell as illustrated in Fig. 2. (The curves of Fig. 6 
show the actual load-deflection characteristics ob- 
tained for the numerical case considered.) 

In Fig. 2, if we consider u as a function of 6, then 
lu, is the value of the indicated maximum at the point 
@ and u is the value of the minimum at point b. The 
portion of the curve oa represents the pre-buckled 
variation of deflection with increasing pressure, while 
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the upper part of the curve, bc, represents the post- 


buckled variation. The portion of the curve labeled 
ab represents a highly unstable condition which cannot 
be realized in the normal physical situation. The value 
“, in Fig. 2 is the upper buckling load and ws is the 
minimum buckling load for the given shell. 

It is now possible to form the envelopes of , and i 
as a function of B and R/t. This relation is given by 


1,2 = Ui »(B, R/t) (47) 


which can be expressed graphically as curves of ™, and 
us versus B for each value of Ft /¢. 

By substituting Eq. (47) into Eqs. (43) and (44), it is 
also possible to graph the deflection and energy at the 
critical values “; and mw as a function of B and R/t. 
This corresponds to the following expressions: 


(6/t)1.2 = D2[B, Ui2(B, R/t), R/t] = Ds(B, R/t) (48) 
Hi. = H.[B, Ui.(B, R/t), R/t] = H,(B, R/t) (49) 


Thus, given a spherical segment with geometrical 
parameters B and K/t, the upper and lower buckling 
loads can be determined by Eq. (47). The deflection 
at the center at the critical buckling loads is given by 
Eq. (48) and the total energy of the system can be ob- 
tained from Eq. (49). 


Discussion of Results— 
The von Karman-Tsien Case 

In the usual discussions (see references | and 6, 
for example) it is observed that upon an increase of 
the external pressure from zero, a clamped spherical 
shell such as shown in Fig. 1 begins to flatten. At 
some critical value of pressure, the shell buckles inward 
causing a deformation at the center in the form of a 
dimple of relatively small solid angle. If the pressure 
is then decreased, a second critical value of pressure is 
reached causing the shell to jump back into approxi- 
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8 ] The curve for g = —1.5421950 is the post-buckled 
solution and g = —0.32584445 is the solution just 
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Fic. 6. Deflection at the center as a function of the load 
parameter. 


mately the form assumed prior to buckling. The pres- 
sure which causes the initial buckling is generally re- 
ferred to as the upper buckling load, and the ‘“‘un- 
buckling’’ pressure is called the lower or minimum 
buckling load. The typical load—deflection curve in 
Fig. 2 discussed previously illustrates this character- 
istic. This figure shows 6, the vertical deflection at the 
center, as a function of pressure for a given clamped 
spherical segment. The pressure 1 defines the upper 
buckling load and wm the minimum buckling load. 
The technique of solution discussed previously yields 
results which clearly reproduce this characteristic. 

We now give the results of the specific numerical case 
considered for the von Karman-Tsien equations. For 
the given value of R/t = 362.5, the solutions, Eq. 
(37), are plotted in Figs. 3a and 3b showing y as a func- 
tion of x. These graphs, which specify the shape of the 
shell, take on various forms depending on the values of 


the parameters w and g. In Fig. 3a, for example, the 


curve for g = —1.575 diverges to minus infinity and the 
second boundary condition is never satisfied. The 
curve for g = —0.S75 takes the form of a damped 


oscillation and, as in the previous case, the second 
boundary condition is not satisfied. The curve for 
q = 0.5 satisfies the second boundary condition at B = 
0.025 and yields an abrupt change in curvature at the 
clamped edge. 

Fig. 3b illustrates the remaining form of solution. 
The four solutions shown in this figure tend to satisfy 
the second boundary condition asymptotically. The 
solid curves are obtained for uv = 0.0011263320, which 
is the urver buckling load. For this value of uw, the 
curve for g = —0.8481995 is the solution just prior to 
buckling and that for g = —1.5658087, just after 
buckling. The dotted curves are obtained for uw = 
0.0010773028 which is the minimum buekling load. 


after “‘unbuckling.’’ The sign of g in the above four 
solutions indicates that the surface of the shell at the 
origin is inverted with respect to its initial shape for 
this specific case. 

For a given uw and R/t, the parameter g determines 
the shape of the shell and thus the value of B [see rela- 
tion (40)]. In the actual machine computation, it was 
found impossible to obtain values of B greater than 
approximately 0.10 (this corresponds to a solid angle 23 
of about 23°). This was due to the fact that arith- 
metic operations on the IBM 650 with automatic 
floating decimal point facility are limited to a precision 
of eight decimal digits. The curves of Fig. 3b were 
obtained only after specifying the parameters ™ and 
qg to one digit in the last place. Computation by double 
precision would undoubtedly improve this situation 
and extend B over a greater range, but this would more 
than double the computing time. Thus, the limitation 
of precision makes it difficult to determine precisely 
the size of the dimple, the deflection at the center, and 
the energy of the system for shells with solid angles 
greater than 23°. 

Fig. 4 is a graph of the relation between & and gq for 
variations in the pressure term u [see Eq. (41)]. — It is 
formed by plotting the values of g which satisfy the 
second boundary condition, Eq. (40), versus the corre- 
sponding values of B. In this figure the curve for 
uw = 0.00120, for example, is single-valued and is typical 
for values of pressure exceeding the upper buckling 
load (uw > m). This curve is asymptotic to g = 4. 
The quantity g, is determined by the condition that 
solutions to the left of g, are divergent and that those 
to the right satisfy the second boundary condition for 
uw = 0.00120. As u is decreased, B as a function of 
q (Eq. +1) becomes multivalued, as indicated for the 
curve « = 0.00113. A further decrease in u produces 
a discontinuity in this function as shown, for example, 
by the curve for « = 0.001125. The gap defines the 
range of solutions, Eq. (37), which are oscillatory, and 
is bounded by two asymptotic solutions. Decreasing 
the pressure w still further causes the gap to widen (see 
curve for u = 0.00109) until « is less than w#. and the 
leftmost portion of the curve vanishes. That is, when 
the value of the pressure is below the minimum buck- 
ling pressure, relation (41) again becomes single- 
valued as is illustrated by the curve for « = 0.00107 in 
Fig. 4. 

Thus, the range of pressure in which buckling can 
occur is determined by the range of solutions which 
cause Eq. (41) to be multivalued. For values of B 2 
0.10, the upper buckling load 1 is found to be the 
value of u which causes the gap to appear (w% = 
0.0011263320). The minimum buckling load i is 
that value of u which causes the curve to the left of the 
gap to vanish (w:, = 0.0010773028). The relation be- 
tween the critical buckling loads and B given by Eq. 
(47) is shown in Fig. 8. This figure is a plot of the 
envelope of 1 and mw for the range of B considered in 
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the computation. The upper and lower buckling 
loads are essentially independent of the solid angle 
of the shell for values of B greater than B = 0.10. For 
0.0575 < B < 0.10, the critical loads increase because 
of the reaction at the clamped edge. The minimum 
value of B at which buckling can occur for the given 
value R/t = 362.5 1s 0.0575, and at this point the upper 


and lower buckling loads coincide, i.e., “4; = ws 
0.00115. 

Eqs. (43) and (44), relating the deflection at the 
center and the total energy of the system, respectively, 
to the quantity 6 for variations in u, are also multi- 
valued in the range of uw in which buckling can occur. 
Fig. 5 is the plot of 6/f as a function of B for various 
values of wu. From this figure we can choose a fixed 
value, B, and obtain values of 6/t for each level of w. 
A plot of 6/t versus u for this B yields the relation 
given by Eq. (45). Fig. 6 is a graphical representation 
of Eq. (45) and indicates the load-deflection char- 
acteristics for various values of B. 

The corresponding graph of Eq. (44) for the energy 
as a function of B for various values of u is not shown 
here because variations of /7 with respect to B are so 
small in the region of interest that the curves tend to 
become confusing even if the scale is exaggerated. 
Proceeding as above, however, we can choose a fixed 
value B and, by numerical interpolation, obtain values 
of IT as a function of u for this B. This leads to Eq. 
(46) which is graphically represented in Fig. 7. 

The portion of the curve for B = 0.080 in this figure 
labeled oa represents the pre-buckled configuration of 
the shell and the portion labeled bc is the post-buckled 
configuration. The part ab represents the highly un- 
stable condition. At the point a, an increment of 
energy A//, is released by the shell as a result of buck- 
ling. Similarly, at the point b, an increment of energy 
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Fic. 8. Envelope of upper and lower buckling loads as a function 
of shell angle opening 


AH, is given up when the shell ‘‘unbuckles.”’ At d, 
the energies of the pre-buckled and post-buckled con- 
figurations are equal but must be distributed differently 
over the shell. The quantity // is negative for the 
range of « shown and represents work done by the pres- 
sure on the shell. 

The relation given by Eq. (47) is the envelope of the 
upper and lower buckling pressures as a function of B. 
This is shown in Fig. 8 which wes discussed previously. 

Relations (48) and (49) corresponding to the de- 
flection and energy at the critical buckling loads as a 
function of B are plotted in Figs. 9 and 10, respectively. 
These functions have been obtained for the range of 
B up to approximately 0.10. The limitation of ma- 
chine precision discussed previously makes it difficult 
to determine these quantities beyond this value of B. 

Eq. (30) of the von Karman-Tsien paper' gives the 
minimum buckling load as 


u = 0.18258 (t/R) (50) 


which corresponds to the quantity m in this paper 
For R/t = 362.5, therefore, the above equation yields 
uw» = 0.00050367 which is about half the value obtained 
by the exact solution (#2 = 0.0010775028). 

Kaplan and Fung® consider solutions in terms of the 


geometrical parameter \ where 
\ = kV R/tsin B (51) 


and k isa constant. 
Convergence difficulties 

values for \ or order 5 or smaller (see page 6 of refer- 

The minimum value of \ at which buckling 


limit their solutions to 


ence 5). 
will occur is given by Kaplan and Fung as 2.1 (see 
Fig. 13 of reference 5). In the present investigation, 
for R/t = 362.5, the maximum velue of 8 obtainable in 
the numerical computation is approximately 11.5° 
which corresponds to a value for \ of about 7. The 
minimum value of 6 at which buckling can occur, 








to 
w 
to 


according to this investigation, is 6.4° which corre- 
sponds to a value for \ of about 4. The discrepancy 
between the minimum values for \, 2.1, and 4, also indi- 
cates the limitation of the power series method. 

Fig. 11 is included in this paper to indicate what 
one would expect from an experimental investigation of 
g as discussed in the introduction. This figure gives 
q, the curvature parameter, as a function of uv, the pres- 
sure term, with a fixed value for R/t. 


Checking the Solution 


One purpose of this investigation is to permit a com- 
parison between the results of the exact differential 
equation given by Eq. (6) and those obtained by the 
von Karman-Tsien' Rayleigh-Ritz method, Eq. (50). 
Since the solutions of Eq. (6) involved certain manipu- 
lations given in the discussion of the method of solu- 
tion, it was desirable to verify that the manipulations 
were performed properly and that the numerical 
quantities were calculated correctly. This was ac- 
complished as follows: Eq. (6) can be rewritten in 
terms of the variables given by Eq. (8). Thus, by 
direct substitution of Eq. (8) into Eq. (6) we obtain 
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Fic. 9. Deflection at the center for the upper and lower buckling 
loads as a function of shell angle opening. 
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where ¢ is a quantity corresponding to the relative 
error of the exact solution of Eq. (6). Now, the com- 
puted quantities x, y, y’, and z’ obtained at each step 
of the numerical integration procedure were substi- 
tuted directly into Eq. (52) and a value was obtained 
for @. This was done automatically by a second ma- 
chine program. The result of this substitution yielded 
an average value of @ on the order of 10~® which can 
be attributed to inherent round-off errors. 

If we had substituted in Eq. (52) values of y’ and 
z’ obtained by graphic differentiation of the y and z 
functions instead of the computed y’ and z’, we would 
have had a test as to how precisely the differential 
equations were satisfied. However, the numerical 
differentiation procedure should not differ significantly 
from the numerical integration procedure and there 
seemed to be no point in testing the integration pro- 
cedure so elaborately. 

However, in the von Karman-Tsien Rayleigh-Ritz 
method, the expressions for y and z are explicitly avail- 
able and can be explicitly differentiated. The above 
check program was applied to these values at each step. 
In this case, the average value of ¢ obtained was of the 
order of 10? or 10° greater than the error in the precise 
numerical solution. This indicates that the Rayleigh- 
Ritz solution given by Eq. (50) is not a solution to Eq. 
(6), the equation of equilibrium. 


Discussion of the Keller-Greenberg-Reiss Case 


A more recent publication by H. B. Keller and E. L. 
Reiss‘ discusses the solution of the system of equations 
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Total energy for the upper and lower buckling loads as a 
function of shell angle opening. 
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given in reference 6 by means of finite differences. 
The step-by-step method of solution described pre- 
viously was applied to these equations in order to com- 
pare this method with the finite difference technique. 
The equations given by Keller and Reiss’ which are 
essentially those given in reference 6 include physical 
assumptions and certain mathematical assumptions 
different from those of von Karman and Tsien.! Their 
equations incorporate tangential and radial displace- 
ments of the shell as proposed by Friedrichs? and 
given by Tsien.* Poisson’s ratio is taken as 0.32 while 
the equations of von Karman and Tsien' assume this 
ratio to be zero. Keller, Greenberg, and Reiss*® assume 
small angles and make other simplifying assumptions 
in formulating the strain displacement relations, and 
these are carried over into reference 7. 

To facilitate comparison with the results discussed in 
the present paper, the equations of Keller, Greenberg, 
and Reiss® were transformed mathematically (but not 
essentially). The geometry of the shell is given in 
terms of the ratio R/t and one half the angle opening, 
6, and the dimensionless pressure term, u, is retained. 
Table 1 relates the symbols used in the present paper 
to those referred to in Keller, Greenberg, and Reiss.® 

The transformed equations corresponding to Eqs. 
(13a) and (13b) of reference 6 are given as 





TABLE 1 


Definition 


This paper Reference 6 





One half the angle opening 


B A 
a 6 The polar angle 
6 a Angle of inclination of the surface of 
the shell 
u a The nondimensional pressure term 
where 
P = 4uR/t Neer 
WW IT Total energy 
p The geometrical parameter where 
P 3(1 — v?) Ro» 
p=2 Jie —B* 
2 t 
¥ A stress function given by Eq. (10) 


of reference 6 





Curvature at the center as a function of the load pa- 
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Fic. 12b. Pre-buckling and post-buckling solutions for a fixed 
shell at various pressures. 
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Fic. 13. Relation between the angle opening of the shell and the 


curvature parameter, for various pressures, 


ay" + 7! = (y/a) + (1/2)(a? — 6) (53a) 


6 R\? ; 
ab” +69’ = (“) + 12(1 — (5) (y8+ ua?) (53b) 
where v is Poisson’s ratio. (Note that the symbols a 
and @ are reversed relative to references 6 and 7.) 
The prime refers to differentiation with respect to a. 
The quantity y is a stress function defined by the fol- 


lowing relations: 
(54a) 
(54b) 


0, = E(y/a) 
eee ee 
o, = Ey 
where o, and a, are the longitudinal and circumferential 
membrane stresses, respectively. 
The boundary conditions equivalent to Eqs. (15) and 
(16) of reference 6 are 


¢6=0 at a=0 (55a) 
6=B68 at a=8 (55b) 
y=0 at a=0 (55e) 
/ Y aaa — a4 
y’—v-=0Oat0=8 (55d) 
a 


Eqs. (53) can be reduced to a system of four first- 
order differential equations by introducing the auxiliary 
variables 

Z = ay’ (56a) 


V = ad’ (56b) 


MARCH, 1961 


This yields the equations 


V' = (6/a) + 12(1 — v?)(R/t)*(y6 + wa?) (57) 


Z' = (y/a) + (1/2)(a? — 6?) (58) 
and 6’ and y’ can be obtained by Eqs. (56). These 


four first-order differential equations are suitable for 
our previous procedures, including the numerical 
integration by the step-by-step method. 
the limit of 6/a as a approaches zero is the parameter 
q corresponding to the curvature of the shell at the 
origin. In addition, the quantities y/a and y’ are 
undetermined at a = 0 and it is necessary to introduce 


As before, 


another parameter, s, such that 


lim y/a = lim y’ = s (59) 

a0 a—0 
The technique of solution is the same as that de- 
scribed previously except that for given values of R/t, 
u, and g, the parameter s must be varied until the 
second boundary condition, Eq. (55d), is satisfied, 
Thus the basic solution can be expressed by the relation 


6 = O(a, g, s, u, R/t) (60) 
which is equivalent to Eq. (37). 

The ratio of the radial deflection to the thickness, 
w/t, as derived from Eq. (11) of reference 6, can be ob- 
tained by integrating the expression 

w'/t = (R/t)(@ — a) (61) 


The dimensionless potential energy as derived from 
Eq. (6) of reference 6 is obtained by integrating the 
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and the angle opening of the shell, for variations in the pres 
sure, 
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6 w 
2v(6’ — 1) (¢ ai 1) — 2| — 4u | (62) 


(Note that the sign of the pressure term, u, is different 
from that in equivalent referenced equation.) 

As an alternative procedure for starting the solution, 
the Runge-Kutta method was substituted for the power 
series method. The solution was continued with the 
step-by-step method previously described. 

The numerical case considered previously was re- 
peated to permit a detailed comparison between the 
finite difference method of solution given by Keller 
and Reiss’ and the step-by-step method. For this 
case, the quantities used were R/t = 362.5 and v = 
0.32. A limited number of computations were carried 
out, the results of which are discussed below. 

The solutions, Eq. (60), which correspond to those 
given for Eq. (37), also take on various forms which 
Fig. 12a, for example, 


describe the shape of a shell. 
indicates that for the given values u = 0.0008 (P =~ 1.4), 
R/t = 362.5, and gq = —1.6, three shapes are obtained 
by considering a as a function of 6, and these depend 
on the value of the parameter s. In this figure, the 
curves for p = 10.51 and p = 19.68 are similar in form 
to those given in Fig. 8b of reference 7. The curve 
for p = 36.75, however, differs in that a slight oscilla- 
tion is observed as a approaches 8. Fig. 12b illustrates 
the slope of the shell as a function of a for p = 10.2 
The form of these curves is also 
There were no such 


for variations in 7. 
similar to Fig. Sb of reference 7. 
divergent solutions as were observed in the von Kar- 
man-Tsien case. This is undoubtedly due to dropping 
various terms in the linearization and simplification 
process. 

Fig. 13 is the equivalent of Fig. 4 of the previous dis- 
cussion. It shows the relationship between g and the 
angle 8 for various values of « when s has been deter- 
mined to satisfy the boundary conditions. When 
g as a function of 8 is multivalued, buckling occurs, 
since for fixed 8 a number of shapes are possible. 

The curves of Fig. 13 can be regarded as level lines 
for a surface which describes u as a function of 8 and gq. 
Such a surface can be used to determine the relation 
between w and gq for a fixed 6 by passing a plane, 8 = 
constant. Notice that this surface has a saddle point 
at approximately 6 = 0.119, g = —0.200, and u = 
0.000956. The lower buckling load is associated with 
the valley in the w surface in the left of the diagram. 

An analogous discussion can be given for the ratio 
of the radial deflection at the center to the thickness 
as a function of 8, as shown in Fig. 14. This figure is 
equivalent to Fig. 5 shown previously. 

Fig. 15 is a plot of the critical buckling loads and is 
in very good agreement with the results of Keller 
and Reiss.’ In the present discussion, for R/t = 
362.5, the relation between p and 8 is given by p = 
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Upper and lower critical buckling loads. (See also 
Fig. 7 of reference 7.) 


Fic. 15. 


841.2436? and the relation between P and uw is P 
1682.49u. Thus, for a given p the value of 6 is deter- 
mined. For this given value of 6 one finds the least 
value of « for which there is more than one value of gq. 
For this least value of u, there will be a positive and a 
negative value of g. The latter value of g is also the 
point where, for fixed u, 8 is a local minimum when 
regarded as a function of g. The locus of such (8,q) 
points, which can be obtained from Fig. 13, determines 
the lower buckling-load curve in Fig. 15. 


In Fig. 13 there is also a locus of points where, for 
fixed 1, 8 is a local maximum as a function of g. This 
locus ties in with a minimum locus at the saddle point 
and determines the curve for the upper buckling load 
in Fig. 15. 

Notice that the accuracy required to establish the 
above phenomena is relatively high. The distinction 
between maximum and minimum points for the uppet 
buckling curve of Fig. 15 is not present in the Keller 
and Reiss’ discussion. 


The results of this investigation yield radial deflection 
characteristics similar to those shown in Figs. Sa and 
9a of reference 7. The relation between the potential 
energy & and the pressure at the critical buckling loads 
is also similar in form to Fig. 6 of reference 7. 


There are considerably more phenomena that could 
be investigated. However, it must be appreciated 
that the simplifications represented by dropping terms 
in the current discussion makes this whole procedure 
The disappearance of the singularities along 
It would 


doubtful. 
with these terms is particularly disturbing. 
be very desirable to repeat this discussion in a mathe- 
matically complete fashion, eliminating, of course, the 
small-angle assumption. 
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Conclusions 
Precise numerical solutions of the von Karman- 


Tsien equation of equilibrium for spherical shells under 
uniform external pressure have been obtained by use 
of the step-by-step method of integration. These 
solutions indicate clearly the configurations of the shell 
before and after buckling, and yield distinct values of 
the upper and lower buckling pressures. These critical 
buckling loads are independent of the angle opening of 
the shell when the angle opening is greater than the 
solid angle of the dimple. A comparison between the 
results of the precise solution and of the Rayleigh-Ritz 
solution indicates that the latter does not satisfy the 
equation of equilibrium. The precise results depend 
critically on terms which are generally neglected in 
the classical linear theory. 


The step-by-step solution process was also applied 
to the mathematical formulation given by Keller, 
Greenberg, and Reiss in order to compare this method 
with the method of finite differences. 
the step-by-step method yielded a finer resolution of 


In this case, 


the buckling phenomena than would otherwise have 
been revealed. The results also indicate the desira- 
bility of a complete mathematical formulation of the 
problem, in which no terms would be dropped and small 


angles would not be assumed. 


Appendix A 
The power series coefficients used in the starting pro- 
cedure are listed below. The coefficients C apply to 
Eqs. (31) and (32) and coefficients K apply to Eq. (33). 
To facilitate the computation, the A’s are given in 


terms of the C’s. 


C; =o 
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Satellite Recovery Techniques for 
Optimization of Touchdown Accuracy 


DANIEL L. ROSAMOND* 
McDonnell Aircraft Corporation 


Summary 


The variables affecting touchdown accuracy of retrorocket- 


recovered satellites of nonlifting configuration are interrelated 
jn ways such that 


For a given nominal orbit there are positions for retro- 


various natural optimum combinations 
exist 
rocket application such that touchdown location is insensitive to 
burnout velocity. There are also special points for minimizing 
effects of errors in burnout flight-path angle and burnout alti 
tude. These optimizing points do not vary rapidly with de- 
partures of the orbit from nominal. 

The natural application of the foregoing is that selecting the 
proper point to apply retrorocket (velocity change) can signi- 
ficantly improve touchdown accuracy. For best overall accuracy 
the chosen point should be based on weighting of the special 
points according to the relative accuracy to which each burnout 
variable can be controlled. Touchdown accuracy is also im- 
proved by choosing the particular retrorocket direction at which 
variations in the direction have no first-order effect. 

Typical data on these effects, obtained by numerical solution 
of the equations of motion, are presented and discussed. A 


brief discussion of the equations of motion is given in an appendix. 


Symbols 
V earth-referenced velocity of mass point representing 
satellite 
V, circular orbit velocity 
climb angle 
h altitude 
AV; velocity change produced by retrorocket 


bp retrorocket direction—angle in orbital plane between 


negative of V just prior to AVr and the vector 


representing AVr 


tr time at which AV» is applied, seconds after burnout 

’ radius from center of earth to mass point representing 
satellite 

¢ latitude, measured from equator 

9 = earth-referenced longitude, measured eastward from 


Greenwich meridian 
¥ azimuth angle of earth-referenced velocity, measured 
from local north 


w angular velocity of earth, 0.0000729211 rad./sec. 

radius of earth at equator, 2.0926 X 10’ ft 

ry radius of earth at poles, 2.0856 & 107 ft. 

fe = radius of earth at latitude ¢ 

G = gravitational constant times mass of earth, 1.4077 X 
10'* ft.*/sec.? 

J oblateness constant, 0.001638, dimensionless 


Other symbols used only for discussion of equations of motion 


are defined in the appendix. 


Subscripts 
BO = burnout 
TD = touchdown 


Presented at the Space Flight Session, IAS National Summer 
Meeting, Los Angeles, June 16-19, 1959. Revised and received 
July 7, 1960. 

* Project Group Engineer. 


Introduction 


Pipe OF THE TOUCHDOWN POINT of recoverable 
artificial satellites to some reasonably small area 
is, of course, an important design consideration, and 
the problems involved in such control have been and 
are being studied quite widely. However, for the most 
part the tendency has been to analyze the accuracy 
requirements and penalties by cascading tolerances 
or errors for separate pieces of the total trajectory from 
launch to 
cussions regarding portions of the total trajectory have 


touchdown. Several analyses and dis- 


been published. Reference 1, for example, divides the 
recovery into three phases: (a) orbit to atmosphere, 
(b) re-entry, and (c) final deceleration; and mentions 
that total range error will be a combination of errors 
in (a) plus errors in (b) and (c). But the graphs and 
discussion presented pertain mainly to effects in each 
phase separately. Reference 2 presents a dispersion 
analysis in which several of the important range error 
contributions are derived and discussed, but neglects 
the effects of the atmosphere even though impact is 
deals solely with the re-entry 


discussed. Reference 3 
phase, and reference 4 deals solely with the phase from 
retrorocket application to re-entry. Reference 5 also 
contains some graphs and discussion of errors and 
accuracy requirements in the retrorocket-to-re-entry 
phase. But the analysis of recoverable satellite touch- 


down accuracy on a launch-to-touchdown basis has 


apparently not been given much attention. This is 
natural, in a sense, because obtaining the desired tra- 
jectory for the launch portion (for example) certainly 
involves problems quite different from those of the re- 
entry. But the piecewise approach is not well suited 
for analysis of touchdown accuracy because it tends 
to obscure or miss altogether the effects of interrela- 
tions of variables in one portion of the overall trajectory 
with those of another portion. 
sider a range of errors of burnout velocity, flight-path 


To be specific, con- 
angle, and altitude. Each of these produces a range 
of errors of re-entry velocity, angle, and location. 
But the combinations of re-entry conditions are not 
arbitrary; e.g., in general the maximum re-entry ve- 
locity does not correspond to the maximum downrange 
location of re-entry. And the correspondence which 
does occur is definitely affected by when and how the 
retrorockets are applied between burnout and re-entry. 
The possibility of choosing a particular retrorocket 
technique which makes at least some of the errors 
produce cancelling effects was recognized from results 
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of studies in which the trajectory from burnout to 
touchdown was analyzed as a whole. It turns out 
that there are several such choices of technique which 
can significantly improve touchdown accuracy. The 
following sections of this paper present data and dis- 
cussion on some of these for a simple mission of a 
satellite recovered after a few cycles of an orbit below 
the radiation belt. There are probably similar effects 
for other types of missions. 


Interrelations Between Burnout, Retrorocket, 
and Re-entry Conditions 


Since the error-cancelling effects are somewhat sur- 
prising and not obvious at first, it is considered worth 
while to include the following qualitative explanations. 


Burnout Velocity Errors 


Consider two orbits which differ only in burnout 
velocity—e.g., those of Fig. 1. They have the same 
velocity at two different points, one on the way up and 
the other on the way down. If the same velocity 
change (retrorocket) is applied to each at the first of 
these points, the high-l’zo9 case will certainly touch 
down farther downrange because it is higher and 
climbing more steeply at the time of the change. But 
for the case of retrorocket applied at the second equal 
velocity point, the effects of the burnout velocity 
difference are compensating—i.e., the steeper dive 
angle tends to offset the higher altitude. In other 
words, orbit differences due to burnout velocity errors 
can produce accumulating or compensating touchdown 
errors depending on where the retrorocket is applied. 
This is shown by Fig. 2, in which the effects are trans- 
lated into differences of re-entry conditions. Note in 
Fig. 2 that for retrorockets applied anywhere within 
about 60 degrees on either side of perigee (assuming 


burnout to be at perigee) the re-entry in the high- 
Vro case is faster, shallower, and farther downrange— 
i.e., all the effects are in the same direction. But for 
retrorockets applied near apogee, the high- go case re- 
enters more steeply. Just after apogee the difference 
in re-entry position falls off sharply, while the steeper 
re-entry effect changes more gradually. This suggests 
that somewhere beween apogee and about 60 degrees 
from perigee the re-entry angle effect may offset the 
re-entry position effect (and also the re-entry velocity 
effect, which is minor). If the retrorocket is applied 
at a selected ¢ime after burnout rather than at a selected 
position, there is also the offsetting effect of the period 
change. The high-Vz 9 orbit has the longer period, 
and hence positions in this orbit lag those of the low- 
Vo orbit. This effect eventually predominates over 
all others because the lag increases with time. To sum- 
marize the effects of burnout velocity errors: If retro- 
rockets are applied shortly after burnout, a higher 
Veo makes touchdown farther downrange; if retro- 
rockets are applied at a certain point after apogee, 
there will be minimum and possibly zero sensitivity 
to Vo errors; and if retrorockets are applied according 
to time rather than position, the lower Vzo makes 
touchdown farther downrange (for a sufficiently long 
time in orbit), and hence there is always at least one 
time for retrorocket application such that touchdown is 
insensitive to Vgo. 
Burnout Flight-Path Angle Errors 

Fig. 3 shows data on two orbits differing only in 
burnout flight-path angle. Note that these orbits have 


common V and / at 180 degrees after burnout and also 
at 360 degrees after burnout. But the climb and dive 
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Fic. 2. Re-entry conditions for Vzo = V- + 150 fps in 
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Fic. 3. Altitude and velocity variations for orbits differing by 1 
degree in yao. 


angles for the orbits at these points are interchanged. 
Hence if the same retrorocket were applied at 180 
degrees, the —ygzo case would touch down farther 
downrange because it is the climbing case at 180 
degrees. Just the opposite would occur if the retro- 
rocket were applied at 360 degrees, and hence at some 
point between 180 degrees and 360 degrees, retrorocket 
application would produce the same touchdown location 
for both orbits. This is the point of insensitivity to 
errors in burnout flight-path angle. One other com- 
pensating effect involving errors in flight-path angle 
is worthy of note: the touchdown errors are not sym- 
metrical about zero ygo. This comes about because 
the perigee altitude is an even function of ygo, but the 
location of the perigee is an odd function of ygo. The 
combined effects thus accumulate for one sign of 
Yro and counteract for the other. Fig. 11 shows this 
nonsymmetrical effect. Note that sensitivity of touch- 
down to ygo errors is, in general, considerably less for 


dives than for climbs. 


Burnout Altitude Errors 

These effects are similar to those of Vg go errors, 
but generally less pronounced. Note in Fig. 4 that 
for orbits differing in burnout altitude there is no veloc- 
ity crossover point corresponding to that for orbits 
differing in burnout velocity. The altitude and climb 
angle differences are also less marked (compare Figs. 
1 and 4). Fig. 5 shows the same general trends of 
re-entry condition differences as shown by Fig. 2, but 
note particularly that the percentage variation of re- 
entry position difference is much smaller for hyo vari- 
ation than for Vgo variation. This means that there 
is less chance of a zero hgo sensitivity position, although 
there will certainly be a minimum sensitivity position 
which also will be located somewhere after apogee 


Da 
/200 
/100 
/000 
900 


800 





ALTITUDE (1000 FT) 








700 
hgo = 820000 FT. 
——-——-— hgo = 780000 FT 


258 ;>——,—__4 
25.6 
254 


25.2 





VELOCITY (1000 FPS) 


250 Eee EE EE EE eee 








O 30 60 90 120 150 180 210 40 270 30 330 %0 


ANGULAR TRAVEL IN ORBIT PLANE 
(DEG. FROM BURNOUT) 


Fic. 4. Altitude and velocity variations for orbits differing by 
40,000 ft. in hgo. 


and before about 60 degrees from perigee (for orbits 
like those used for illustration). For retrorocket firing 
according to time rather than position, there is the addi- 
tional effect of period change due to burnout altitude 
difference similar to that due to burnout velocity dif- 
ference. And likewise, this effect eventually pre- 
dominates so that there will always be at least one time, 
although not necessarily on the first orbital cycle, where 
touchdown is insensitive to burnout altitude. 
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Fic. 5. Re-entry conditions for hgo = 
of increments from conditions for hgo = 780,000 ft. AVr = 
500 fps. 
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Retrorocket Angle Errors 

There is always a particular value of 6,, the retro- 
rocket direction angle, in the first quadrant such that 
touchdown is insensitive to 6g variations about this 
point. Here again the combination of an even func- 
tion effect with an odd function effect produces the 
cancellation of errors. The tangential component of 
the retrorocket velocity change is an even function of 
dz while the normal component is an odd function. 
Hence the effects of these components on touchdown 
will combine for one sign of 6g and oppose for the other. 
Plus values of 5g give combining effects, i.e., downward 
and backward thrust increments both tend to move 
touchdown farther uprange. But in the first quadrant 
(+6,) the variations of the components due to vari- 
ations in 6g have opposite signs; increasing dz de- 
creases the tangential component but increases the 
normal component. If both components affected the 
touchdown equally, the optimum would be 45 degrees, 
where the gradients of the components with 6d, are 
equal and opposite. However, the tangential com- 
ponent usually has the greater effect, and this causes 
the optimum 4d, to be less than 45 degrees, where the 
tangential gradient is less. 

The explanation of error-cancelling effects at opti- 
mum 6g can also be made on the basis of variations of 
the re-entry conditions. Fig. 6 shows that as dz is 
varied from 0 to 10 degrees there is practically no 
change of re-entry angle or re-entry airspeed (which is 
because the tangential component is practically un- 
changed), but there is a marked uprange shift of the 
re-entry position (which is due to the normal com- 
Thus, from zero dz to 10 
But from 40 


ponent at 10 degrees 6p). 
degrees 6g touchdown shifts uprange. 
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veral. Burnout: V = V. + 100 fps; y = 0; & = 800,000 ft.; 


Cp (S/m) = 0.6. 


degrees dp to 50 degrees 6g there is practically no 
change in re-entry position (because the effects of the 
component variations are about equal and opposite). 
Hence in this range the effects on re-entry angle and 
re-entry speed will predominate. A higher 6, makes 
the re-entry faster and shallower (because of the re- 
duced tangential component), and hence from 40 de- 
grees dz to 50 degrees 6g touchdown shifts downrange. 
This is opposite to the effect from 0 to 10 degrees dp, 
and hence at some point between these two ranges, 
touchdown will have zero gradient with dp. 


Discussion of Typical Data 


For detailed analysis of touchdown position vari- 
ations, including specific data on all the effects dis- 
cussed in the preceding section, an extensive family of 
numerical solutions of the differential equations of 
satellite motion, including effects of earth rotation and 
oblateness, has been obtained by use of an IBM 704 
these equations is given 


computer. Discussion of 


in the appendix to this paper. 


Nominal Orbit 

Fig. 7 shows the path, in terms of latitude versus 
longitude, of an orbit corresponding to a straight-east 
launch from Cape Canaveral. This orbit is simply a 
typical one selected as nominal for the purpose of com- 
puting touchdown variations due to errors from the 
nominal. Burnout is nautical 
miles from Cape Canaveral, and the nominal burnout 
conditions are Vgg = V,. + 100 fps, yao = O, and 
hgo = 800,000 ft. All of the data on touchdown posi- 
tion are presented in terms of longitude only, because 
all the touchdown points computed fall within a narrow 
band near the path of the nominal orbit. Errors which 
cause an orbit period change do, of course, shift the path 
because of the rotation of the earth, but this is a small 
effect for touchdowns after only a few orbital cycles. 
Even after a day in orbit (17 cycles) the maximum 
path shift due to Vgo errors is only about one nautical 
mile per foot per second, and this shift occurs only 
near the equator; it falls off to zero near the maximum 
latitude of the orbit. 
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Optimum Retrorocket Angle 

Fig. 8 shows the effect of 6g variation on touchdown 
longitudes for a AV’z of 500 fps applied at three times 
covering a wide range of orbital positions. The same 
general effect is shown for all three fg values, but there 
is a slight shift of the optimum 4, point. This shows 
up better in Fig. 9, where the slope is plotted versus 
5g for fg = 3,000 seconds and tg = 5,000 seconds. 
The curve for tg = 4,000 seconds falls so near the one 
for 5,000 seconds that it is omitted for clarity. Actu- 
ally, Fig. 9 shows practically the whole range of vari- 
ation of the slopes versus 6, because 3,000 seconds after 
burnout is not far from apogee and 5,000 seconds is not 
far from perigee. Hence Fig. 9 shows that for AVz = 
500 fps the optimum 6, for any point on the selected 
nominal orbit is confined to a range of less than 10 de- 
grees, roughly between about 30 degrees and 40 de- 
grees. There is also a variation of optimum 6, with 
the magnitude of the velocity change, as shown by 
Fig. 10. A larger Al’z makes the optimum 6, higher. 
For the rest of the IBM computations, i.e., for investi- 
gation of effects of variables other than retrorocket 
angle, 6g was fixed at 40 degrees and Al’, was fixed 


at 500 fps. 


Effect of Burnout Errors—First Pass Touchdown 


Fig. 11 shows a summary of touchdown locations for 
fairly wide ranges of burnout errors, each varied singly. 
The y range of plus or minus one degree may seem small, 
but in terms of the vertical velocity component it is 
over 400 fps rate of climb or dive. The insensitive 
points are readily apparent from these plots. Note 
that except for some of the y variations, the gradients 
of touchdown position with burnout errors are reason- 


ably constant. In other words the slopes, and hence 
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errors. AVr = 500 fps, Cn (S/m) 0.6. 

the points of greatest error-compensating effect, do 
not vary rapidly with departures of the orbit from 
nominal. Also, these points are not particularly criti- 
cal, i.e., errors are small for a considerable region near 
the optimum. Of course, there really is some ‘“‘cross 
coupling”’ between the variables—e.g., the point of zero 
y slope does depend somewhat on Igo and hgo, and 
similar dependencies exist for the other optimum 
points. But these effects are fairly small. A 100 
fps change of Vy shifts the time for zero sensitivity 
to y errors only about 60 seconds. 
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nominal burnout: V = V, + 100 fps; y = 0; hk = 800,000 ft.; 
Cpo(S/m) = 0.6. 


Note in Fig. 11 that applying retrorockets at 3,000 
seconds after burnout (slightly after nominal time for 
apogee) is unfavorable for accurate touchdown because 
of the steep slopes for all three burnout variables. If 
the amount of retrorocket is reduced so that touchdown 
for tg = 3,000 sec. is at the same longitude as touch- 
down with AV, = 500 fps at a more nearly optimum 
tp—say 4,500 sec.—then the near-apogee case is even 
worse because the effect of AV, errors is greater at 
lower AVz. A detailed study, reference 6, of just 
such a comparison of apogee firing versus near-opti- 
mum-time firing showed that the optimized technique 
provides about a 3 to 1 reduction of the 99 per cent prob- 
able touchdown area. 


Cyclic Variation of Sensitivities 

It has been pointed out that there are times or posi- 
tions for retrorocket application such that natural 
compensation of touchdown errors due to burnout 
errors results, and that there are other times or posi- 
tions causing accumulation of error, sometimes up- 
range and sometimes downrange. Actually this is 
a consequence of the cyclic nature of the variation of 
the partial derivatives of touchdown position with 
respect to the burnout variables. Fig. 12 shows these 
cyclic variations for the first three orbital cycles. 
These results were obtained by numerical differenti- 
ation of touchdown data computed for variations 
each way from the nominal burnout conditions. 

From this summary one can see at a glance which 
times are favorable for applying retrorocket and which 
are unfavorable. Of course the overall best and worst 
times depend on the launch guidance, i.e., on the rela- 
tive accuracy to which each burnout variable can be 
controlled. Perhaps there are optimizing schemes 
which could be exploited during launch also so that 
burnout errors would occur only in combinations in 
which the sign of one of the partial derivatives is 
opposite to that of the other two. But from Fig. 12 
it is apparent that such combinations, and hence the 
launch guidance scheme used to produce them, would 
be different for different retrorocket application times. 

Note also from Fig. 12 that the V and / derivatives 
show a one-way shift as well as a periodic variation. 


AEROSPACE SCIENCES 


MARCH, 1961 

This is due to the period change effect discussed in a 
previous section. The y derivative does not show a 
one-way shift because y variations do not change the 
Because of the one-way shift due to 
applying retro- 


orbit period. 
period changes, the technique of 
rockets according to a preselected time after burnout 
becomes unsuitable after about 6 or 8 orbital cycles. 
Beyond this the period change effect becomes so large 
that there are no zero sensitivity times for Vgo and 
hgo errors, and still later even the least sensitive times 
have prohibitively high partial derivatives with respect 
to Vgo and hgo. Fig. 13 shows these partial deriva- 
tives for a portion of the seventeenth cycle. Note that 
here the least sensitive time is worse than the most 
sensitive time for the first cycle. Thus for recovery 
after a large number of orbital cycles, the retrorocket 
must be applied according to position rather than time 
so that the period change effect will be eliminated. 
However there are still favorable points and unfavor- 
able points—i.e., the cyclic effects still cause optimums, 
even though the optimum positions may not be zero 


sensitivity positions. 


Appendix: Equations of Motion of a Point Mass 
Moving with Respect to a Rotating Earth 

The equations of motion which were integrated 
numerically to obtain the results presented in this 
paper are as follows: 


¥ = re’? + (rcos? ¢)(6+ ww)? +4A,4+ g, 


¢ = —27%~/r — (cos gsin ¢)(6+ w)? + 
(1/r)(A, 2 £y) 
6 = 2(¢gtan ¢ — #/r)(6 + w) + (1/r cos ¢)Ag 


where A,, A,, and A, are the radial, northward, and 
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eastward components of aerodynamic force per unit 
mass, and g, and g, are similar gravitational com- 
ponents. There is no g, component for a rotationally 
symmetric earth. The rest of the terms come from 
expressions for acceleration components in spherical 
polar coordinates. Derivations of these expressions 
can be found in textbooks, e.g., reference 7. 

Resolution of the 4 terms from the usual L and D 
(lift and drag) can be done in terms of a bank angle B, 
which locates the L vector in the plane perpendicular 
to the velocity V, and the azimuth and climb angles, 
y and y, which orient the V vector relative to the local 
radial-northward-eastward frame. It follows from 


simple projections that 
A, = (1/m)[L cos B cos y — D sin y] 
A, (1/m)[—L sin B sin y — 
(L cos B sin y + D cos y) cos yw} 


Ag (1/m)[Z sin B cos y — 
(L cos B sin y + D cos y) sin y] 


where m is the mass of the moving body. The com- 
ponents of V are #, rg, and r6 cos ¢ in the radial, north- 
ward, and eastward directions, respectively, and hence 


V = VP + re? + £6 cos? ¢ 
y = arc sin (#/V) 
y = arc tan (6 cos ¢/¢) 


Lift and drag are computed from the usual nondi- 
mensional coefficients C; and Cp based on the reference 
area S and the dynamic pressure '/2:pV*. For this 
computer program the air density p was taken to be 
a function of altitude only, with the altitude calculated 


from 
h=r-r,=r1-—tfre V cos? go + (re/rp)* sin? ¢ 


For altitudes below 280,000 ft. the 1956 ARDC model 
atmosphere was used, but for the higher altitudes p was 


—5.74 
56356) | 


which gives better agreement with existing satellite 


obtained from 
p = 1.8036 X 107% X 


1 + h/20860000 


orbit decay data. 
The gravitational terms in the equations of motion are 
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orbital cycle AVr 500 fps, dr 40°. Cp (S/m) 0.6. 
determined by differentiation of the potential function 
for an oblate ellipsoid, which for small oblateness is 


approximated by 


U = (G/r){1 + (J/3)(re/r)?11 — 3 sin? ¢g) +... | 
g, = OU /dr = -—(G r) [1 + I re r)? x 
(1 — 3 sin? ¢)] 
gp = (1/r)(QU/O¢g) = —(G/r*)J(re/r)? sin 2 ¢ 
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Recent Studies on the Effect of Cooling on Boundary- at the end of the potential core of the jet. The nozzle used in — 
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RICHARD J. WISNIEWSKI AND JOHN R. Jack 250 the supercritical case was designed to produce a shockless jet. 
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The investigators were familiar with Sprenger’s ideas on resonance 
tubes and suggested that the measured temperature rise might 
result from unsteady heating processes within the stagnation 
chamber of the total-temperature probe. 

A third aerodynamic energy separation which can be brought 
to mind is that of the Ranque- Hilsch vortex tube, in which cooled 
ir is produced at the core and warmed air along the walls of the 
In practical cases the flow in the chamber 
s highly turbulent. The most plausible explanation of the 
so-called Ranque effect is that it is due to turbulent mixing 
wross the large radial pressure gradients within the vortex 
chamber, the radial motions of the expanding and contracting 
as essentially adiabatic and as being 
This concept 


yortex chamber. 


lumps being conceived 
followed by mixing with the surrounding fluid. 
was first put forward in connection with the vortex tube by 
Knoernschild.? The calculations of outlet temperatures made 
by van Deemter* and by Deissler and Perlmutter® using theo- 
retical models based on this idea gave results which agreed well 
with measured values. <A similar hypothesis has proved useful 
in problems of turbulent mixing in the atmosphere.’ 

Sprenger tested a Hartmann generator and a Ranque tube of 
the same pressure ratio, finding their 
On this 


similar dimensions at 
energy-separation performances to be very much alike 
ground he advanced the theory that the Ranque phenomenon 
was due to unsteady effects like those he held responsible for the 
tubes. this was due 


But the opposite conclusion—that the 


heating of resonance Apparently idea 
originally to Ackeret 
resonance-tube separation mechanism can be explained in terms 
f that of the vortex tube—seems equally convincing. 

In view of the evidence available it seems more plausible to 
isctibe both of these effects (and that noted by Glassman and 
John® in free jets) primarily to turbulent mixing, and to suggest 
that the unsteady processes within the resonance tube may well 
be of secondary importance. Turbulence seems to be the one 
factor common to all these energy separations; shock systems 


and organized unsteady flows are not universally present. 
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Sloshing of Liquids in Spherical Tanks 


James D. Riley and Nat. W. Trembath 
Space Technology Laboratories, Inc., El Segundo, Calif. 
July 29, 1960 


 ypegeeed PAPER by Budiansky! discussed sloshing of liquids 
in spherical tanks. A program? to compute the sloshing in 
axially symmetric tanks of quite general shape is available at 
Space Technology Laboratories, Inc. 

This program applies a Rayleigh-Ritz procedure using three 
deep-tank eigenfunctions and five odd powers of the radius 
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fffvorav, J foras 
V ‘ Ss 


The triple integral is extended over the volume, V, occupied by 
the fluid at rest; the double integral is extended over the quies- 
cent free surface, S. 

The program computes the three lowest eigenvalues and the 
corresponding wave shapes together with certain constants de 
scribing a mass-spring analogy. The program will handle ring 
shaped tanks and tanks with maximum tank radius greater than 
the radius at the free surface 

The program has been run for a spherical tank at various fluid 
levels with the results shown in Table 1. These results provide 
curves (see Fig. 1) agreeing quite well with the conjectured curves 
of Budiansky,' and also the experimental results of McCarty and 
Stephens.5 

The amount of time on the IBM 709 for this program was ap- 
proximately two minutes setup time for a particular tank shape, 
and one minute to determine the output for three modes at each 
surface level considered. 

In theory, being a Rayleigh-Ritz procedure, the program gives 
upper bounds for the eigenvalues. Of course, numerical evalua- 
tion of the integrals involved invalidates the theory, due to round- 
off. However, a Rayleigh-Ritz procedure developed by Henrici® 
using spherical harmonics leads to rational coefficients for the half 


TABLE | 
Eigenvalues for a Spherical Tank. 
(The notation is that of reference 1.) 





é Ni Ne As 
0.8 4.32 9.81 15.43 
0.6 2.77 7.06 11.14 
0.4 2.14 6.00 9.53 
0.2 1.79 5.50 8.79 
0.0 1.56 §.27 8.49 

—0.2 1.39 5.25 8.57 
—0.4 1.26 5.40 8.99 
—0.6 1.16 5.73 9.97 
—0.8 1.07 6.29 12.18 
—0.9 1.03 6.65 14.33 
—0.95 1.02 6.84 

—0.98 1.01 6.94 

—0.99 1.003 

—0.995 1.002 

—0.999 1.0003 
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sphere. This procedure gives 1.5618 in the 4 by 4 case, and in v = Ov/dt = 0, t 0 (3d 
the 6 by 6 case a value of 1.5602 as upper bounds for \;._ In the : 

general program, the value obtained to four decimal places was We suppose that 

1.5594. This indicates that Budiansky’s value of 1.565 was ac- v= 7, + 2% { 


curate to within 0.3 per cent. 
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Note on Forced Torsional Oscillation of a 
Circular Cylinder by the Application of a 
Transient Force at One End 


A. K. Mitra* 
Department of Mathematics, Jadavpur University, Calcutta, India 
August 1, 1960 


ee OBJECT of this note is to solve the problem of torsional 
vibration of a circular cylinder of which one end is fixed 
while the other end is acted on by a transient shearing force, 
the curved surface being free from stress at all times. It is 
believed that the problem has not been attempted by any previous 
investigator. 


SOLUTION OF THE PROBLEM 


We use cylindrical coordinates 7, 6, z, the z-axis coinciding with 
the axis of the cylinder. 

We suppose that the end z = 0 of the cylinder is fixed and the 
end z = 7 is acted on by a transient shearing force proportional 
to re, the curved surface being free from stress at all times. 
We suppose uw; = uz = 0 and v(= ug) is independent of 6. 

Then two equations of motion are identically satisfied, and 
the remaining one is reduced to! 

p Iv O77 4 1 ov 4 ra) (1) 
pw Ol? or? f OF ra) ; 
Also we have 


- 


rr = 0, 996 = 0, 53 = 0, 75 = 0 


We are to obtain a solution of (1) subject to the conditions 


g=0, s=0, #2>0 (3a) 


, 


E v al 
‘ia | ¥ibews ; 


me = Sre™, 14>0 (3c) 
HI =I 


(S being a constant), and 


t>0 (3b) 





* The author wishes to express his sincere thanks to Professor B. Sen, 
D.Se., F.N.1., of Jadavpur University for his kind guidance in the prepara- 
tion of this note. 


where v7 is a particular solution of Eq. (1) satisfying the condi 
tions (3a), (3b), and (3c) while vp is the function that satisfies 


(1) and the conditions (3a) and (3b), and also we suppose that 


M[Ov2/Oz]2 = 1 0 


Let 
: of - 
%, = gr sinh kse“ (5 
Now (5) satisfies (1) if 
b2 _ (p/p)Q? 02 /¢2 (6 
where 
-2 - 
c? = p/p (7) 
Again, as (5) satisfies (8c) we have 
g = S/uk cosh kl (8 
We assume 
te =f p> Sin Gn2(Cn COS ant + Dy sin ant) (9) 
n = 0 
As (9) satisfies (1) we have 
° 2. 2 
Qn° = C°Gn* (10) 
Also from the relation w[Ov./0z]z .7; = 0 we have 
2n+1 7 
g 11) 
») Ml 
Again we have 
iv] 20 = z + v»] 2-0 = 0 
from which it follows that 
=C, sin 9.8 = —g sinh ks 


l l 
therefore C, f, sin? g,2dz = -¢ f, sinh kz sin q,2dz 


S8elk cosh k/ 


therefore C, = (-1)"* 
$k2]/2 + (2n + 1)?x? 


Also we have 


H 7 [= 4 =| i 
| a t=0 ot ot _|i= 


from which it follows that 
l . l . . 
anDn . sin? g,2dz = ¢Q " sinh kz sin g,2dz 
Therefore 


ott - ink 8/Q¢k cosh kl 7 (13) 
4]2k? + (2n + 1)?22 


So we get ultimately 


= 


¢ 


aoe sinh kze~@ + 
uk cosh kl - 


8rSl = (—1)” 


wn o mide? + (2n + 1)2x?} 


sin (2m + 1)rz sin ant , 
—————— | = O08 at +} (14 
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Flow Induced by a Cavity in a Supersonic 
Stream 


H. Serbin 

Senior Staff Engineer, Hughes Aircraft Company, 
Culver City, Calif. 

July 21, 1960 


a ERRATIC BEHAVIOR of a free or partially supported body 
in the presence of an open bomb bay is probably due to the 
aerodynamic disturbances generated by the bay. In a practical 
case, the open bay displays various kinds of stores and equip- 
ment so that the shape is very irregular. However, for a first 
step, it is convenient to idealize the bay as a long, rectangular 
cutout in an infinite plane (Fig. 1). 

For the sake of simplicity, the calculations will be carried out 
on the basis of linearized aerodynamics.!:2 The underlying 
assumption that the flow deflections created by the body are 
small cannot be true either at the front or rear of the bay. These 
are local areas, however, and it is possible that the large features 
of the flow might be predicted, at least qualitatively, by the 
linearized theory. 

The perturbation potential g(x, y, z) due toa body, symmetrical 
relative to the plane y = 0,* is 


7) S gy(&, 0, $) x 


y,s) = -(1 


dedt/V(x — ¢)? — B(x — £)? + y%], B? = M?—1 (1) 


Let y = J(x,s) represent the bounding surface near the 
plane of symmetry. Then 


g,(x, 0, 3) = U(OY/dz) 
where U is the velocity in the free stream. For a cavity** of 
depth h, width 2a, and the length c 
Y = —hS(s) + hS(z — c) “| << 
= 0 x| >a 


The rear edge of the cavity will have no influence on the flow 
in front of that edge and will, therefore, be disregarded in the 


following. Then* 
g(x, O, 2) = —hU8(z) Ixl <a 
= 0 lx] >a 
Eq. (1) becomes with this substitution 
g(x, y, 3) = (hU/r) I dt V3? — B[(x — £)? + y?] 


where the last integration is to be carried out overall, — on the 


interval —a < & < a for which the points (£, O, O) lie within 
the retrograde Mach cone issuing from the point (x, y, 3). In- 


& 


tegrating on & yields 


g(x, y, 3) = (hU/7B) : sin~! (B tan a) (2 
a 
where tana = (~ — x) IV 22 — B*y?, ltana! < 1/B (3) 


and the summation is carried out over the a (or £) corresponding 
to the end points of the £-interval lying within the cone. 

It is convenient to consider separately the regions shown in 
Fig. 2. Region 0 consists of points (x, y, z) whose retrograde 
cones do not contain any points £. For such points ¢ vanishes 
Region I consists of points whose cones are inter- 
Then, from Eq. (2) 


identically. 
sected twice by the é-interval. 
o(x, y, 3) = (hU/xB)[(2/2) — (—2/2)] = hU/B 

* The subscript y represents partial differentiation. 


** The Dirac delta-function below is represented by 65(x), the unit step 


function by S(x). 
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Fic. 3. Surfaces of constant potential 


In both regions, one can use the combined representation 


Region0 + J: g(x, y,2) = (hU/B)S(z — By) (4) 


The perturbation velocity is therefore*** 

Vo = —hAU|j — (1/B)k]-6(z — By) (5) 

Points of regions II and II’ (Fig. 2) have the property that 

the retrograde cone cuts the £-axis only once. Considering only 
region II, Eqs. (2) and (3) give 

+ 2/2)] (6) 


g(x, y,2) = (hU/xB)[sin= (B tan a 


where tan a = (a — x)/V2? — By? (7) 


Within this region, the perturbation velocity is 


hU 
Veo = - > 4 
™V 22? — B2[(x — a)? + y?| 
,. Bia — x)y, (a—x)z .] 
=i + v es ke (8) 
| 3? — Bry? 2? — Bry? f 


The potential function, Eq. (6), depends only upon a, hence 
the velocity vector is normal to the surfaces tan a = constant. 
The intersections of these surfaces (looking upstream) with a 
plane normal to the free-stream vector are shown in Fig. 3. 
The arrows show the direction of flow 

It readily can be seen from Egs. (5), and (8), and the usual 
formula for perturbation pressure that the pressure induced by 
the bay is zero within region I, positive in that part of region II 
above the bay, and negative in the rest of region II. 





*** The presence of the 6-function is due to the approximation implied 
by the thin-surface theory, according to which all singularities on the 


surface of the body are moved to the plane of symmetry 
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Fic. 4. Schematic of streamlines in elevation. 


The flow in region III is obtainable by superposition of the 
flow in II with the antisymmetric flow in II’. The details 
will not be discussed further. 

On the basis of the theory above, it is clear how air flows into 
the cavity. As an air particle moves over the leading edge of 
the cavity, it passes through the 6-function zone (actually an 
expansion wave followed by a shock wave). The particle is 
suddenly deflected downward and then is set into motion along 
the bay in a manner depending on the position of the particle 
relative to the Mach The shown 
schematically in Fig. 4. 


cone. various cases are 
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Comments on ‘Stability of Flat, Simply 
Supported Corrugated-Core Sandwich Plates 
Under Combined Loads’’ 


Paul Seide 

Technical Staff, Space Technology Laboratories, Inc., 
Los Angeles, Calif. 

July 12, 1960 


point out a misleading 
Auelmann 


HE PURPOSE OF THIS NOTE is to 
concept in reference 1. Messrs. 
are under the impression that they have made the assumption 
that any line in the sandwich core that is initially straight and 
normal to the middle surface will remain straight after de- 
formation of the plate but will deviate from the direction of the 
normal to the deformed middle surface by an amount that is 
proportional to the slope of the plate surface, this proportionality 
factor being constant throughout the plate. However, the 
method of solution presented in their paper appears to be equiva- 
lent to solving for the shear angle in terms of the deflection 
function, so that, within the limits of approximation of a 
truncated Rayleigh-Ritz approach, their solution is an exact 
solution of the Libove-Batdorf sandwich equations. 
Let us consider the third of Eqs. (2) of reference 2, with the 
longitudinal shear stiffness De, of the core assumed infinite. 
Then in the notation of reference 1, 


' L—-p.D oO D — p? a oy 
£ toe U oy? Cy = Oy \Ox? oy? “ 
(1) 


Harris and 


If w is chosen in the form 


wi) 


— - _ mae 8 
w= 5 ) Qmn Sin sin ¥ (2) 
a b 


m=l1n=1 


and Q, in the form 


ub 
m=ln=1 


, max 
sin ~ COS 
a 


nry 


b (3) 


mn 


we have, from Eq. (1) of the present paper, 
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L (nx/b)|(mx2/a)? + (nx/b)?]D 
Imn 


1+ [(1 — p)/2](D/U)(mx/a)? + (D/U)\(nr b)? 


If we define 


OW mn OW m, Owe 
Lmn . — 5 
O7 oy l 
we have 
Om l l Riis 
mn 1 — i-_ = 
OW mn/ OV U (nz/b)am, 
2/)(a?/b?) — (1 + pw)m? 
= (6 
2n*(a?/b?) + (1 — p)m? + 2J(a?/b?) 


which is identical with Eq. (10) of reference 1. Thus, if the 
authors have used a different value of g in each equation of their 
stability criterion, corresponding to different values of m and n 
{this is not clearly discussed in the paper), they have, in effect, 
assumed Q, to be given in the form of Eq. (3) and have mini- 
mized the potential-energy expression with respect to the coefi- 
cients b,, to obtain an exact solution to the problem, rather 
than the approximate solution they believe they have. 

For the case of a simply supported plate under longitudinal 
compression, with the core longitudinal, the results of references 
1 and 2 are identical. The differences indicated by Table 3 of 
reference 1 may be due to the different values of Poisson's ratio 
used in the two papers (0.3 in reference 1 and 1/3 in reference 2 
and to inaccuracies in plotting and reading the charts of reference 
2. 

It is interesting to note that a method such as was originally 
contemplated by the authors was used in reference 3. Since an 
approximate deflection function with only one undetermined 
was assumed, the determination of the constant ¢ 
The use of the method with deflection 


coefficient 
was relatively simple. 
functions having more than one undetermined coefficient was 
not indicated in reference 3 and it is not at all certain how one 
would proceed analytically to minimize the buckling coefficient 
with respect to the constant g and obtain a workable solution 
One can, of course, assume various values of g and minimize 


graphically. 
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Note on the Sufficient Conditions for the 
Stability of Axially Symmetric Flows to Axially 
Symmetric Disturbances 


G. N. V. Rao* 

Department of Aeronautical Engineering, 

Indian Institute of Science, Bangalore, India 

August 1, 1960 

em LINEARIZED PERTURBATION DIFFERENTIAL EQUATION 
for the problem under consideration is 


d? 1d gy” 7 
(U — C) ee —ar®lot+ — UYU" og = 
dr? r dr r 


* It is a pleasure for the author to acknowledge the interest and encourage 
ment that he has received from Professor Dhawan during the course of this 


work. 
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jin the notation of reference 1. The purpose of this note is to 
bring out some of the similar properties of axially symmetric 
and plane flows, as far as sufficient conditions for stability are 
concerned. As an example we shall consider the flow over a 


slender body of revolution in axial motion. Only the 


very 
incompressible case is considered. Taking 7 (7? — 1)/2, 
Eq. (1) can be transformed to 
’ . i d? 4 reer 
(l C)i (1 + 29) — — a? Io — (1 + 2y) U6 = 
dn? 
asi ” sedi , a 

(1 + 2n)? > + 4(1 + 2n)? ee ——- a ap () 

iaR U] dn dn? 


Multiplying Eq. (2) by ¢/(1 + 2”), where ¢ is the complex 
conjugate of ¢, and integrating between 0 and o, we get 


[2 + 2a?h? + ato? —1aRQ + taRC(I? + a*ly?) (3) 


where 
= $6 . “—* 
Iy? =. an, I? = ¢'¢'dn, 
0 (1 + 2n) 0 
@ o''o"" 
J.2 = <a dn (4) 
g (1 + 2n) 
and 
“a I , , a? rer ‘ rr , 
Q = U¢'d’ + >, + (U" #bo Idn + U'o'odn 
J0 1 + 2n 0 
(5) 


From the solution of the differential Eq. (2), it is known 
(to be published by the author; also, reference 5), that it behaves 
as e~” for large values of ». Hence following Lessen? one con- 
cludes that the integrals in Eqs. (4) and (5) are bounded. Adding 


to Eq. (3) its complex conjugate, one obtains 
iaR(Q — Q) + 2aRCiW(T? + aly?) = 
—2(Io? + 2a7t;? + ato?) (6) 


where Q is the complex conjugate of Q. 

This equation is similar to Eq. (3.3.2) of reference 1, obtained 
for the case of two-dimensional parallel flows. Hence all the 
analyses of references 2 and 3 can be extended to flows with axial 
symmetry with little modification. It also shows that in view 
of the fact that the solutions of Eq. (2) behave similarly to the 
well-known solutions of the Orr-Sommerfeld equation, it is reason 
able to expect similarity in the pattern of the stability loops 
of the two cases, where they both exist. This statement, how- 
ever, has to be qualified by the following considerations. 

The disturbances considered are axially symmetric. Exami- 
nation of the linearized equations of motion for the disturbances 
indicate that an azimuthal disturbance which is independent 
of the azimuthal coordinate (of the type encountered in the 
well-known Taylor-Gortler vortices) behaves independently of 
the longitudinal and lateral disturbances, and in such a cir- 
cumstance an analysis based only on the axisymmetric dis- 
turbances is valid. Reference 4 in addition also shows that 
such azimuthal disturbances are always stable disturbances. 

The necessity for using only axisymmetric disturbances is 
dictated by the extremely complicated equations of motion for 
a general three-dimensional disturbance. There is also the 
difficulty with regard to Squire’s theorem which has served so 
well in two-dimensional parallel flows. Examination of the 
complete equations of motion, following Squire’s approach, 
seems to indicate that this theorem is not valid for the case 
of three-dimensional disturbances in axisymmetric mean flows. 
The author has learned from a personal communication that 
Mr. Schade of Berlin has also reached similar conclusions. 

Similarity of the general nature of sufficient conditions for 
stability does not mean that a flow in the axisymmetric case 
behaves similarly to its counterpart in plane parallel flow, for 
which the well-known example of the plane and axisymmetric 
Poiseuille flows can be given; but the former view has been 


FORUM 249 


advanced by observation of the stability loops calculated by 
Tatsumi® for the flow in the inlet length of a circular pipe 

For examples other than the one considered here it is not 
difficult to arrive at the proper normalizing coordinate to bring 
Eq. (1) into a form similar to Eq. (2). For example, 7» = 
(1 — y?)/2 would be suitable for pipe flows and r?/2 would 


be suitable for jets, to lead to equations similar to (2) 
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Sub-Alfvenic Flow in Magnetoaerodynamics 
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T° REFERENCE 1, McCune and Resler considered steady, plane, 
magnetoaerodynamic flows past thin cylindrical bodies with 
uniform, parallel magnetic field in the undisturbed part of the 
stream. Within this category they treated a subcase of particular 
interest, namely, ‘‘aligned-fields’’ flow, in which the undisturbed 
magnetic field is parallel or antiparallel to the free stream. 
For this subcase there are some unique phenomena; for ex- 
ample, if the electrical conductivity is infinite the velocity and 
magnetic vectors are everywhere parallel and the field vanishes 
inside the body, assuming only that it is not also a perfect con- 
ductor. If the conductivity is large but finite, there is a bound- 
ary layer of large current density and vorticity, outside of which 
the inviscid-perfect-conductor approximation holds. The mag- 
netic-field strength inside the body is then not exactly zero, but 
is small. 

The purpose of the present note is to discuss the exterior (in- 
viscid, perfect-conductor) flow; namely, to propose certain ex- 
planations of McCune’s and Resler’s results and, incidentally, to 


correct some errors in their paper 


THE TANIUTI-RESLER DIAGRAM 
It is shown in references 1 and 2 that the equation describing 
small-perturbation, plane, aligned-field flow of a perfect con- 
ductor is 
(1 — M*)\ver + yy = 0 (1) 


where y is the perturbation stream function, such that -—y, = 


v’ and yy, = (1 — M..?)u’, u’ and v’ being the perturbation ve- 
locity components, 7. the Mach number of the free stream, and 
M is defined by 
Mm? = M.2m?/(M..2 + m? — 1) (2) 
where m is the Alfvén number, i.e., the ratio of the stream speed 
U to the Alfvén-wave propagation speed. (This ratio is called 
1/A, in reference 1.) If H. is the magnetic-field strength in the 
free stream and p, the density there, the Alfvén-wave speed is 
H../V4xp.. Thus, M is a kind of effective Mach number, and 
U/M is the speed of propagation, in the free-stream direction, 
of the acousti-mhd waves that replace both sound and Alfvén 
waves in a compressible perfect conductor. 
Fig. 1, which also appears in reference 1, shows the areas of an 
m, M., plane in which Eq. (1) is hyperbolic and elliptic, respec- 
tively. The same diagram was used by Taniuti* in a study of the 
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Fic. 1 (/eft). Diagram showing elliptic and hyperbolic (cross- 
hatched) regimes of steady plane flow of a perfect compressible 
conductor with aligned fields. Fic. 2 (right). Hyperbolic sub- 
sonic plane flow of a compressible conductor past a flat-plate air- 
foil. 


characteristics of aligned-field flows. It pertains not only to 
small-perturbation flow but also, locally, to plane flow without 
linearization. 

Let us discuss briefly the three regions, A, B, and C, in which 
m <1, i.e., sub-Alfvénic flow. 
Region A: m? + M.,.? < 1. Here MU? < 0. 
simple affine transformation of the Prandtl-Glauert type, Eq. 
(1) is transformed into Laplace’s equation. The only unusual 


By means of a 


feature here is that the Prandtl-Glauert factor V1 — M? is 
greater than 1, so that the distortion of coordinates and body 
profile is opposite to that to which we are accustomed. There are, 
however, no obvious qualitative changes compared to conven- 
tional gasdynamics: pressure variations will be qualitatively 
conventional, ete. 
Region B: Here the flow is hyperbolic, but, as pointed out in 
reference 1, the ‘‘Mach waves” that make it up are inclined up- 
stream. Thus the pattern of streamlines and waves is as shown 
in Fig. 2. Now, this is subsonic flow; hence, speed reduction is 
indicated by wide spacing of streamlines, and speed increases by 
narrow spacing. Thus, for example in Fig. 2, so far as the pres- 
sures outside the current layer are concerned, positive lift and 
drag will occur on airfoils, ete. 
Region C: Here M < 1, and although the flow is supersonic, it 
is elliptic. It is related to irrotational incompressible flow by a 
conventional Prandtl-Glauert transformation. But this only 
means that the streamline pattern is affinely related to incom- 
pressible flow; one must not forget that (1 — J/..2)u’ = dy/dy, 
and (1 — M.,?) < lhere. (Although this relation was correctly 
written in Eq. (13) of reference 1, Eq. (23) of that paper is in- 
correct.) Thus, the speed variation with stream-tube area is 
typically supersonic in this régime: speed is greatest where 
stream-tubes are widest, and vice versa 

We, therefore, have the interesting combination of a typically 
elliptic streamline pattern with supersonic speed-pressure vari- 
ations. This implies that the pressure is least near stagnation 
points and greatest where the related incompressible flow has 
lowest pressure. Unless surface pressures are reversed in sign 
by current-layer forces, we must expect negative lift at posi- 


tive incidence, etc 


NET PRESSURE AT BODY SURFACE 

To ascertain the net pressure at the body surface, and thence 
lift and drag (if any), we must correct the pressures discussed 
so far for the body forces on the current layers. Qualitatively, 
this is very simple, for the streamlines and magnetic lines are 
coincident in these flows. Thus the streamline pattern tells the 
story of the field intensity. We need only adopt the rule: 
Where stream-tubes are narrow, field intensity is great, and vice 
versa. This rule, together with Ampére’s Law, suffices to de- 
termine the sign of the net surface currents—e.g., into the paper 
in Fig. 2. A moment’s reflection will confirm that this is the 
sign of the net current (upper surface minus lower) for a lifting 
airfoil in all the cases discussed so far—viz., in Regions A and C 
because the streamlines are carried over from incompressible 
flow, and in Region B because the streamlines are as shown in 
Fig. 2. Thus the lift due to surface currents, its sign being 
given by j X H, is negative in all three cases, in disagreement 
with the statement below Fig. 2 of reference 1. 

Alternatively, the rule can be stated even more simply by use 
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TABLE 1 





‘Lift Due Lift ] ue 
to to 


Streamline External Surface Net 
Region Spacing Flow Pressures Currents Lift 
A Conventional Subsonic + - - 
subsonic 
B As in Fig. 2 Subsonic -+ _ + 
_& Conventional Supersonic — - = 
subsonic 
D - Subsonic - _ + 
E Conventional Supersonic + + - 
supersonic 
of the concept of magnetic pressure. The streamline spacing 


tells us the magnitude of this pressure, whose expression is H2/87. 
If the magnetic pressure is greater above than below the airfoil, 
as in Fig. 2, a negative contribution to lift results. 

These qualitative results are collected in Table 1 above. 
For completeness, we include Regions D and E, using the same 
two rules-of-thumb—viz., (a) magnetic lines are spaced like 
streamlines—hence lift-due-to-surface-currents is negative when 
stream-tubes are narrow above airfoil and wide below; (b) pres- 
sure variation within stream-tubes is determined solely by 
whether flow is subsonic or supersonic. Finally, we also include 
the sign of the net, or total lift for positive incidence. This 
naturally depends (except in C and E) on the relative magni- 
tudes of the two contributions, and is always given correctly 


by Eq. (22) of reference 1. 
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Recent Studies on the Effect of Cooling 
on Boundary-Layer Transition at Mach 4 


Richard J. Wisniewski and John R. Jack 

Lewis Research Center, National Aeronautics and 
Space Administration, Cleveland, Ohio 

August 2, 1960 


pew ADVENT of high-speed flight has necessitated the study of 
boundary-layer transition on highly cooled bodies. In- 
vestigations such as those of references 1—4+ have concentrated 
on this problem and have indicated, contrary to the trends 
predicted by small-disturbance theory, that premature transition 
can be found with cooling. This phenomenon, commonly called 
‘transition reversal,’’ has been investigated and is discussed in 
detail in references 2-5. 

The purpose of this note is to report some recent transition 
data obtained on a cooled cone in a Mach 4 wind tunnel. The 
model, a sharp-tip cone (included angle 13.5°), was cooled by 
liquid nitrogen to a temperature of —340°F. The cooling 
method and the data analysis are similar to that described in 
reference 3. 

The data are analyzed in terms of the ratio of wall temperature 
to adiabatic wall temperature and the local transition Reynolds 
number based on the surface distance to transition. Shown in 
Fig. 1 are the present data and the flight results of reference 6. 
The wind-tunnel data at a Reynolds number per foot, Re/ft. 
of 7 X 10%, indicate the expected variation of the transition 
Reynolds number Re;,, with cooling. These data show no 
“transition reversal’? even though temperature ratios as low as 
Tw/Taw ~*~ 0.22 were investigated (minimum value for these 
tests). In fact, for values less than 7,/Ta» ~ 0.52, the entire 
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body is covered with laminar flow. These data may be inter- § Jack, John R., Wisniewski, Richard J., and Diaconis, N. S., Effects of 
a: preted as indicating a region of complete stabilization below Extreme COE OF B seen Layer Transition, NACA TN 4094, 1957 
. / 0.52. tl } : ll with the limi f 1 4 Stetson, Kenneth F., Boundary-Layer Transition on Blunt Bodies With 
= (.! -reby agreeing well w > . te “e : : ; 
A I IZ, neve V agreeing We we 1 the ~ ot! complete Highly Cooled Boundary Layers, Journal of the Aero/Space Sciences, Vol. 27 
Net stabilization predicted by the theory of reference 7. No. 2, p. 81, February, 1960 


Jack, John R., and Wisniewski, Richard J., The Effect of Extreme Cooling 


ift However, when the Re/ft. is raised to a value of 9 K 108 the 
and Local Conditions on Boundary-Layer Transition, Journal of the 


i. preceding clear-cut effect of cooling is not found. Rather, we : é, 
Aero/Space Sciences, Readers’ Forum, Vol. 25, No. 9, pp. 592-593, Sep 


















find that Re,, increases until a value of T,,/Ta, = 0.54 is obtained, Pi egg 
+ and below this temperature ratio a peculiar effect of cooling is 6 Rabb, Leonard, and Disher, John H., Boundary-Layer Transition at 
i“ noted. As the temperature ratio was decreased from approxi- High Reynolds Numbers as Obtained in Flight of a 20° Cone-Cylinder With 
iad 5 f. ll c ocal : ré € perc R¢ os ¢a le *f 255 5 
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a ° ~ os cca ~ é = | u . aicuiattons ré lé itty of the aminar oundary 
of Retr is normally referred to as “transition reversal. How- Layer in a Compressible Fluid on a Flat Plate with Heat Transfer, Journal of 
— ever, as the temperature ratio was decreased further to approxi the Aeronautical Sciences, Vol. 19, No. 12, pp. 801-812, 828, December, 1952 
mately 0.35, a totally unexpected result was found. The * Wisniewski, R. J., A Comparison of Boundary-Layer Transition Data 
, a — 1 F 7 Ra sis ° if : 1 ti . Obtained in Free Flight With Data from Laboratory Sources; to be published 
‘Ing transition teynolds a — reversed itsell — eaenar in 1960 Proceedings of the Bumblebee Aerodynamics Panel, University of 
Sr. flow was established over the entire body. This double reversal Michiran. Ann Arbor, May 9-10, 1000 
foil, and the region of complete stabilization, found below 7y/Toaw ~ 
0.35, represent a new phenomenon not reported previously. It ; ‘ 
ve. should also be noted that again (as in, e.g., reference 5) the Re/ft. 
ume as well as the temperature ratio appears to play an important Solutions to the Heat-Conduction Equation 
like role in the transition picture. With Time-Dependent Boundary Conditions 
hen Although the fairing of the curve connecting the three vertical 
res- points may appear to be arbitrary, it is actually based on a 
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Fic. 1. Effect of cooling on boundary-layer transition. ** Professor of Mechanical Engineering. 











252 JOURNAL OF THE AEROSPACE SCIENCES—MARCH, 1961 
i 
f { . 
1s 
| 
(N,X) = DS BCX** Zoi 4 541 (12) 


————— 
Q(X) 








Fic. 1. One-dimensional heat conduction in a slab 


027(N, X)/ON?2 — O7r(N, X)/OX + K +MN? = 0 (2) 


The initial temperature distribution is assumed to be represented 
by an even-order polynomial of the form 


7(N,0O) = FN? + GN (3) 
The boundary conditions are as follows: 
07r(0, X)/ON = 0 (4) 
O7(1, X)/ON = Q(X) (5) 
where the heat flux in Eq. (5) can be expressed as a polynomial 
of (s + 1) terms, 
QX) = > H,X* (6) 
s=0 


Eqs. (2)-(6) were solved by means of the Laplace trans- 
formation. The solution for the temperature distribution is 
given by 


8 
V(N,X) = SO WX*Zoay — (2F + 4G)Z, — 24GXZ; + 
s=0 
(2F + K)X + (12G + M)X* + (M + 12G)XN? + 
FN? + GN* (7) 


where 
© 
Zosy1 = 2*t1g1X' SO [i+ erfe [(2r + 1 — N)/2X"“] + 
r=0 


i+! erfc [((27 + 1+ N)/2 t/a) } (8) 


In solving aerodynamic-heating problems, it is often con- 
venient to specify the separate components of the flux, the 
conductance ratio, and the temperature potential. Assume 
the conductance ratio can be represented as a polynomial of 
(4 + 1) terms, 


CxX)= >} cx' (9) 
i=0 


and assume that the temperature potential may be represented by 
ro X) = (taw — tw) = BX (10) 
where 7 is an even integer. The heat flux then becomes 
; 
Q(X) = DF BCX‘ (11) 
i=0 


If the initial temperature distribution is uniform and if the heat 
generation vanishes, the solution for the temperature distribution 


7=0 


The transient temperature distributions represented by the 
simple polynomials of Eqs. (7) and (12) can be readily computed 
with the aid of the design charts.!. These design charts present 
values of Z(X, N) for a wide range of the parameters found in 
actual practice. The general formulas and the design charts 
make it possible to calculate rapidly transient temperature 
distributions for such cases as: (1) aerodynamic heating of a 
thin wing, (2) stagnation-point heating of a blunt nose cone, 
(3) cooling of an unclad fuel element in a nuclear reactor, (4) 
heat treating and quenching of slabs, and (5) mass-transfer 
problems such as drying. 
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Some Simplifications in the Treatment of 
Rotary-Inertia Effects for Transverse 
Vibration of Beams 
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A“ EXTENSION of the method of influence coefficients to 
include the effects of rotary inertia in the transverse 
vibration modes of beams was first introduced by Scanlan.'? 
We describe briefly this treatment of the problem, in which the 
deflections y and slopes y’ are defined through the integral 
equations. 


b b 
Wx) = of f br(x, £)m ()y() dé +f. d(x, e)ile)y' eae | 
a a 


(1) 
i ; b b ; , 
y'(x) = ol ff or(x, &)m(E)y(E)dE + f. om(x, &)i(E)y (=)ae | 
a a 


where w is the natural frequency, m(x) the running mass of the 
beam, 7(x) the running moment of inertia per unit span, about 
the local chordwise axis; and the integrals extend over the 
length of the beam. From the static theory of beam deflections 
the results summarized in Table 1 are readily obtained. 
It can also be seen that 
or(x, £) = (0/dx)or(x, £) | 
(6) 


om(x, £) = (0/0x)bm(E, x) ] 


bu(x. &) = or(é, x) 
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The equations (1) can be treated numerically by approximating 
the integrals with a suitable quadrature formula or by lumping 
the masses at properly chosen points. In either case one is led 
to a matrix-eigenvalue problem in which the order of the matrix 
is double the number of points chosen to approximate each 
integral. It is however possible to derive from (1) a single 
integral equation which involves only y’ and consequently 
requires for its solution the treatment of a matrix whose order 
is the same as the number of masses 

Since Eqs. (1) hold in the strict sense only for cantilever-type 
motion, we consider this case first. Taking the equation for 


y’, namely 


b b 
yix) = | f or(x, £)m(£)y(=)dé +f oul x, £)i(€)y"(E)dE | 
a a a 
{ 


} 


ind integrating the first term on the right by parts gives 


i b b b 
y'(x) = w*, y(a) or(x, E)dé 4 y'(E)dé or(x, 2)m(2)dz 
Me a a g 
b 
+ om(x, £)y’"(E)i(E)dE (8) 
a 


For a cantilever beam y(a) = O and the first term on the right- 
hand side vanishes. The result is a single integral equation 


b 
+ f. or(x, z)m(s as | yg ae | 
é 


for y’ 


g b> 
y(x) = ot Dl | om v, E)i(E 
a 


(9) 
For a free-free beam the equations (1) must be rewritten: 
eb 
W(x) = yo + xbo + w? J bp(x, £)m(E)y(E)dE + 
a 
b 
bu(x, &)1(E)y"(E)dE 
a 
b 
y(x) = do + w? or(x, &)m(E)y(E)dE + (10) 
a 


b 
f. om(x, £)i(E)y(E ae | 
a 


where yp and ¢» define the rigid-body motion of the beam. These 
unknown constants yo and @» are determined by the conditions 
that the sums of all forces and of all moments must be zero: 


b 
m(x)y(x)dx = 0 
a 
b b “/ , ‘ 
m(x)y(x)xdx + i(x)y'(x)dx = 0 
a a 


If the beam has symmetry, it is convenient to measure x from 


(11) 


the center and to consider separately the cases where the de- 
flection is either symmetric or antisymmetric. In the first case 
% = 0 and the moment equation is satisfied a priori, while in 
the latter case yo) = O and the force equation is satisfied auto- 
matically. In the general case it is convenient to take x as 


measured from the center of gravity, so that 


b 
m(x)xdx = 0 
a 


The constant yo is then found by multiplying the first of Eq. (10) 


by m(x) and integrating, giving 


b b 
—-My = of m(£)y(£)dE f br(x, £)m(x)dx + 
a a 
eb b 
J i(é)y'(é)dE Su(x, £)m(x)dx (12) 
a a 


b 
M = f m(x)dx 
a 


Similarly, multiplying the first equation by [m(x)-x] and the 


where 


9 
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second by i(x), integrating, and adding gives 


\ b 
—I¢ = of m(£)y(t)dE X 


b b 
6r(x, £)m(x)xdx + op(x, &)i(x)dx + 
a a 
eh b 
i(é)y(E)dE bu(x, E)m(x)xdx + 
a a 
eh ) 
om(x, &)i(x)dx ( (13) 


b *h 
J = m(x)x*dx + 1( x )dx 
wa a 


These values when substituted into Eqs. (10) lead to a new 


pair of equations for the free-free modes: 


b b 
(x) = wif P(x, &)m(£)y(E)dE 4 ff. Q(x, Epil e)y' (a | 
a a 
b eh 
y'(x) = w? P(x, &)m(é)y(E)dé + f Q.(x, &)i(E)y"(E)dE 
a a 


(14) 


with 


As before, the second of these may be converted into a single 


equation for y’: 


- \ b 
y'(x) = w? ) y(a) PAx, 2)m(s)dz 4+ 
a 


or “bh ) 
Q.(x, &)u(E) + J, P(x, 2)m(2)dz y(E)dE¢ (15) 
a . 


The constant y(a@) may be found by making use of the condition 


for equilibrium of forces: 


b 
f y(~E)m(E)dE = O (16) 
a 


Upon integrating by parts, this condition becomes 
b s 
My(a) + M(é)y(E)\dE = O (17) 
a 
where 
b 
M(x) = f m(&)dé, M = M(a) 
z 


Thus the kernel of the integral equation for y’ becomes 


b 


Q.(x, &)i(E) + P(x, z)m(2z)dz — 


Me) f° 
P(x, 2)m(2)dz (18) 
M 


a 
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Concerning Real Gas Effects on Hypersonic 
Inlet Performance 


C. W. Hartsell 

Member, Technical Staff, Space Technology Laboratories, Inc., 
Los Angeles, Calif. 

August 17, 1960 


SYMBOLS 


Cp = specific heat at constant pressure 
h = enthalpy, altitude 
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Fic. 1. Inlet path indicated on an air Mollier diagram. 
KE kinetic energy 
M Mach number 
P pressure 
R gas constant 
S entropy 
T temperature 
U velocity 
n kinetic energy efficiency 


Y ratio of specific heats 


Subscripts 


r = total or stagnation 
2 = station at end of compression process 
3 = station at end of the isentropic expansion from station (2) to the 


ambient pressure P , 
« = free stream 
PG = perfect gas 
RG = 


s = isentropic process 


real gas 


a FUTURE of air-breathing propulsion for hypersonic flight 
speeds depends to a substantial extent on the attainable 
The present supersonic-flow analytical and 
part 
owing to the relatively unknown influence of the real gas effects 


inlet performance. 


experimental techniques are of questionable value, in 


on inlets. In order to shed some light on this problem area, 
the simplest inlet (i.e., normal shock inlet) will be evaluated by 
both real- and perfect-gas techniques and the resulting per- 
formance compared. 

An elementary use of the first and second laws of thermo- 
dynamics shows that the increase in entropy during an adiabatic 
compression process may be written as follows: 
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Fic. 2. Hypersonic normal-shock inlet performance comparison 


(real gas vs. perfect gas). 
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AS = Cp lin (T?7,/Tr.) — Riu (Pr./Pr..) (perfect gas) (1 


For a perfect gas it is known that 77, T7..; thus the entropy 


rise is solely a function of the inlet pressure recovery. However, 
if we consider real gas flow where Cp, R, etc., are dependent on 
T and ?, it is obvious that the entropy rise is some function of 
two state variables such as JT and P. Thus in comparing inlets, 


pressure recovery alone is an inadequate parameter. The 
fundamental concept of kinetic energy efficiency (7) may be 
used to describe the inlet performance at hypersonic flight 
speeds. 7 is defined as the ratio of the kinetic energy of the 
captured air (KE) when it is isentropically re-expanded to the 
ambient pressure divided by the free-stream kinetic energy 


(EZ). 


computations 


Eq. (2) presents the formulation used for real-gas 


KE, _ U;? hr. —h 


RG KE, U2 (hr. — hs)s 


To appreciate nzg more clearly, note Fig. 1 where a typical inlet 
process is indicated. 

With the assumption of a perfect gas, Eq. (2) may be manipu 
lated into the standard form as given in Eq. (3): 


[Pr./Pr.| 
npg = 1 - = ss 1.40 3 


The use of reference 1 along with the Mollier diagram of reference 
2 provides sufficient information to evaluate npg with the use of 
Eq. (2). 
combined with Eq. (3) yield npg 

The results of the above computations are presented in Fig. 2 
A study of the 
namely, (a) the lower the 


The usual perfect-gas normal shock equations when 


as the ratio nrg/npq versus M,, and altitude. 


figure indicates two basic trends, 
altitude of flight the higher the ratio, and (b) as VW 
It should be noted that npg is not dependent on 


increases 
the ratio drops. 
altitude. 

The altitude sensitivity noted comes about due to the effect 
of pressure on dissociation. As the altitude drops the process 
pressure level increases, tending to suppress dissociation, and 
thus the process tends to approach the perfect-gas results. The 
trend noted with increasing J/,, indicates that the real-gas 
detrimental effect on inlet performance becomes more important 
as M 
formance tends to the perfect-gas solutions. The limitations 
indicated 


increases. Conversely, as M/,, decreases the inlet per- 


of reference 2 restrict the low M,, computation as 
From other results in reference 3, it is estimated that the ratio 
tends to unity near V,, = 6 at 100,000 feet. It should be 
noted that the normal shock inlet investigated herein probably 
yields a maximum of real-gas effects. A higher efficiency inlet 
would tend to reduce the real-gas effects due to the higher 
pressures developed which tend to suppress dissociation and thus 
minimize the errors due to assuming a perfect gas. Reference 3 
provides a value of the ratio for an optimum (from shock struc 
ture considerations) two-shock inlet of 0.80 at AZ, 26.7 and 
h = 100,000 feet. This value is shown in Fig. 2 to indicate the 
substantial reduction of the real-gas effects as predicted pre 
viously. 

The results of this investigation lead to two conclusions, 
namely, markedly lower efficiencies due to real-gas effects and 
a first-order altitude effect on efficiency. Both these 
predicted trends are detrimental to air-breathing propulsion 
The first effect indicates that in order to obtain the 
perfect-gas performance level a more efficient inlet must be used 
Secondly, in order to maximize the inlet performance a low 
altitude is preferable, which then tends to maximize the aero- 
As noted previously, a higher efficiency inlet 


inlet 


systems. 


dynamic heating. 
will tend to minimize the above trends. 
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3 Hartsell, C. W., A Study Real Gas Inlets at Hypersonic 
Speeds; submitted as a term paper to the Mechanical Engineering Depart 
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University of Southern California, June 


A Note on Jacobi’s Method for Real Symmetric 
Matricest+ 


Shan S. Kuo 
Assistant Professor of Civil Engineering, Yale University, 
New Haven, Connecticut 
August 12, 1960 
I’ MANY PROBLEMS dealing with the eigenvalues and eigen- 
vectors, it is convenient to use a well-known method due to 
Jacobi! for diagonalizing real symmetric matrices. It consists 
of performing a sequence of plane rotations (orthogonal trans 
formations), each of which reduces to zero one of the nondiagonal 
elements of the matrix, until the sum of squares of nondiagonal 
elements fails to decrease. Less attention has been paid, how- 
ever, to the total number of rotations required to diagonalize 
matrices of various orders 
Recently, in connection with a study of the effect of the 
the 
approximating the dynamic characteristics of a uniform beam,? 


selection of number of masses and their distribution in 


Jacobi’s method was used to calculate the natural frequencies, 


modes, and other related matrices for a range from 2 to 20 


lumped masses using two different methods of lumping (cases a 
through d) as shown in Fig. 1. 


1 The numerical work was done at the Massachusetts Institute of Tech 


nology Computation Center 
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gem o-6—6—4 (a) 
| HJ H H H 
N 49 3 2 i 


Fic. 1. Types 
of lumping and 
support. 
4 3 2 | 
———a—— i) 
H 


EPL By 


B.S 
Ow 
nm 


H Lif 4 fH 


10 20 30 40 60 80100 


ae 


200 300 500 800 
ere Vy 





TTT TY 


20 
15 


TTTTTTTITNTT 


Matrix 


pe ee eee eres 
° 


OrpvER oF 
| 
py 


it 





ereriae iri awe nes ree eT 
30 40 60 80100 200 300 500 800 


fe Seeeeerer oe! 


10 20 








Averace NuMBER OF ROTATIONS 


Relation and of 


between order of matrix number 
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Fic. 2 
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FORUM 2: 


TABLE 1 


Number of Rotations in Jacobi’s Method 


With a series of symmetrical matrices whose elements are 
either the influence coefficients A the 


as in cases d, b, and c, or 


La. | 
matrix 








elements formed by the matrix | Vm Vv mM | as in case d 
(where mj; the Table 1 
the number of rotations for convergence to eight decimal places 
The numerical work was done on an IBM 704 computer with a 


program written in FORTRAN floating-point system 


are elements of mass gives 


In Fig. 2 is shown the relation between the order of matrices 
and the average number of rotations 

Numerical work also indicates that for both cases of hinged- 
l(b) the normal-mode 


functions are the ordinary trigonometric sines; 


hinged beams shown in Figs. l(a) and 


hence they are 


devoid of any discretization error due to the basic approximation 


of replacing the distributed-mass problem by the discrete model 
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The Lienard-Chipart Criteria for the Stability 
of Polynomials 


William Squire 
Senior Research Engineer, Department of Applied Mechanics, 
Southwest Research Institute, San Antonio, Tex. 


September 7, 1960 


r | MHE PROBLEM of determining whether a polynomial is stable 
whether all the roots have negative real parts, is im 
The Routh-Hurwitz 


i.e., 
portant in several branches of engineering 
stability criteria are well known and discussed in many standard 
texts. However, as can be seen from a recent note by Saunders! 
and the references cited there, some important improvements due 
to Lienard and Chipart? are not widely known and are being re- 
discovered An discussion, 
with some extensions of the work of Lienard and Chipart, is 


by various investigators extended 
given in a recent book by Gantmacher* which should be con 
sulted for details which cannot be given in a brief note 

The basic theorem is that the condition for 


f(s = Za, 2" = V0 


to be stable can be reduced to any of the following four sets 
being positive: 
(1) ao, the even coefficients and the even-order Hurwitz deter- 


minants; 
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(2) ao, the even coefficients and the odd-order Hurwitz deter- 


minants; 
(3) ao, the odd coefficients and the odd-order Hurwitz deter- 


minants; 
(4) a, the odd coefficients and the even-order Hurwitz deter- 


minants. 

Therefore, the number of determinantal conditions are cut in 
half. Since the usual criteria of all coefficients and all Hurwitz 
determinants being positive are necessary, satisfying any of the 
reduced set of conditions implies the positivity of the rest. 

It is also possible to compute the number of unstable roots 
from expressions which only utilize the odd or even Hurwitz deter- 
minants, but these are too lengthy to reproduce here. 

The advantages of the simplified criteria in analytical studies 
are obvious, but the decrease in numerical computation in a case 
where the coefficients are given numerically might be questioned. 
One cannot avoid computing the (m-1)th-order Hurwitz deter- 
minant (from which the xth order follows by a single multipli- 
cation). In a systematic scheme, such as that presented by 
Duncan,’ based on expansion in minors, all the lower order 
determinants are evaluated in the course of the calculation. To 
derive any advantage from the elimination of the even (or odd) 
determinants, some other scheme for the numerical evaluation 
of the determinants must be used. Such a scheme might be 
based on the use of eliminants as discussed by Frazer and Dun- 
can,® but further investigation is required to decide on the best 
formulation and to determine if there is an appreciable saving 


in numerical work. 
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Expressions for Damping Matrices in Linear 
Vibration Problems 


P. Lancaster 
Lecturer, Dept. of Mathematics, Univ. of Malaya in Singapore 


August 15, 1960 


— RECENT PAPER by R. D. Milne! on this subject prompts 
me to make this contribution. Milne makes the assumption 
that the damping forces in the problems investigated can be 
represented in terms of the ‘“‘complex’”’ damping concept. Inci- 
dentally, a much more appropriate and descriptive terminology 
for this phenomenon is hysteretic damping (discussed by Neu- 
mark,? the author,* and others). 

The assumption that damping forces in an aircraft structure 
are predominantly hysteretic in type would be very difficult to 
justify. Viscous and nonlinear effects will almost always be 
present as well. However, we can go one step further and as- 
sume that both hysteretic and viscous damping forces are present. 
In theory, at least, we can obtain expressions for the two corre- 
sponding matrices using essentially the same data as Milne. 
That is, if the system has m degrees of freedom (to a semirigid 
approximation) we require 7 complex vectors representing sinusoi- 
dal applied forces, p; say; at given frequencies w,; together with 
the 2 corresponding complex vectors defining the resultant steady- 
state responses, q, Say. 

If A and C are real (” X nm) matrices of inertia and stiffness 
coefficients, respectively, and B and G are real (n X n) matrices 
of viscous and hysteretic damping coefficients, respectively, then 
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the ‘‘response equations”’ are 

(—Aw,? + tBu, + iG + C)qr = pr 
for r = 1,2,...,n. If we define complex (n X n) matrices P and 
Q to have columns p),ps, Pn and qi,Q»,. . .,Qn, respectively, and 
if 2 is the real (n X n) diagonal matrix with diagonal elements 
@1,@0, ,wn, then the above equations can all be rolled into one: 


—AQO? + iBQ2 + iGO + CQ = P 


If we now write P Pr + iP;and Q Qr + iQ1, where Pp, 
P7,Qr and Q; are real, then on equating real and imaginary parts 
of the left- and right-hand sides of the above equation we obtain 
BQ;2 + GQ; = Cr — AQRQ? _ Pr\ 

> ( 1) 

and BOr® + GQr = —CQ; + AQ? + P,S 
The first of these with B = 0 is essentially the result used by 
Milne. However, assuming that Qe and Q; are nonsingular the 


above equations can be solved for B and G to give: 


B= (DrQ 14 DiQr 1) Q72Q0) 1 — QCrOQQVR a 
G= (DrQ 10) 14 D;Q Or VQ 10; 1 _— Or 1Or 1)-1 


where D Dr + iD; = CQ — AQQ? — P 


I feel that any results of this type can only be of academic in- 
terest if applied to evaluating matrices B and/or G, either of 
whose elements are small compared to those of A and C. There 
is some virtue in Milne’s example, in which the effect of omitting 
some degrees of freedom is investigated theoretically, but the pres- 
ent state of the art of resonance testing coupled with the compli- 
cated structures often tested mean that in practice Q cannot be 
measured to very great accuracy. On the right-hand side of 
Eq. (1) therefore, we must subtract two inaccurate large quanti- 
ties in order to obtain something of much smaller order! 

From the theoretical point of view then, it may be of interest 
to state some similar results for systems with viscous damping 
only. These results give A~', B, and C™! (if it exists) in terms of 
the vectors giving the natural modes of free vibration of the 
damped system. For simplicity, and in the spirit of this context, 
we assume that A, B, and C are symmetric and that a// the zeros 
(A,) of 

| Ad? + BA+ C| 


arise in complex conjugate pairs. We form an (” X m) diagonal 
matrix A from these, taking one root from each pair. We define 
T to be the (n X mn) complex matrix whose columns are the vec- 
tors t, defining the free modes of vibration corresponding to the 
selected values of \. That is, we have 

(Ad,? + Br, + C)t, = 0 


for y = 1,2,...,n. It is proved in reference 4 (theorem 3) that 


if the vectors t, are normalized so that 


t,'(2Ar, + B)t, = 1 


then 
A-'= TAT'+TaT’ and CC! = —TA™T’ — TaT' 
where 7 and A are the complex conjugates of T and A. The 
methods of that paper can be used to obtain the damping matrix 
in the form 
B = —A(TA®T’ + Ta2T’)A 


In practice it will be even more difficult to measure the vectors 
t, than the response vectors q,. However, the theory is there 
should the experimental techniques ever become available. 
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at frequent intervals. The work of editorial preparation will be 
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MatTTER Usuatty DELETED: Photographs or illustrations of 
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be clearly designated by name in the margin of the manuscript. 
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